Under consideration for publication in the SIAM Journal of Applied Dynamical Systems 



1 



The Stability and Dynamics of Localized Spot Patterns in 
the Two-Dimensional Gray-Scott Model 

W. CHEN, M. J. WARD 

Wan Chen; Department of Mathematics, University of British Columbia, Vancouver, British Columbia, V6T 1Z2, Canada, 

(Currenntly at OCCAM, Oxford University, Oxford, U.K.) 
Michael J. Ward; Department of Mathematics, University of British Columbia, Vancouver, British Columbia, V6T 1Z2, Canada 

(corresponding author) 

{Received 10 September 2010) 

The dynamics and stability of multi-spot patterns to the Gray-Scott (GS) reaction-diflusion model in a two-dimensional 
domain is studied in the singularly perturbed limit of small diffusivity e of one of the two solution components. A hybrid 
asymptotic-numerical approach based on combining the method of matched asymptotic expansions with the detailed 
numerical study of certain eigenvalue problems is used to predict the dynamical behavior and instability mechanisms of 
multi-spot quasi-equilibrium patterns for the GS model in the limit e — s- 0. For e — >■ 0, a quasi-equilibrium fc-spot pattern 
is constructed by representing each localized spot as a logarithmic singularity of unknown strength Sj for j — 1, . . . ,k 
at unknown spot locations Xj £ Q, for j = 1,. . . ,k. A formal asymptotic analysis is then used to derive a differential 
algebraic ODE system for the collective coordinates Sj and Xj for j — 1, . . . , k, which characterizes the slow dynamics of 
a spot pattern. Instabilities of the multi-spot pattern due to the three distinct mechanisms of spot self-replication, spot 
oscillation, and spot annihilation, are studied by first deriving certain associated eigenvalue problems by using singular 
perturbation techniques. From a numerical computation of the spectrum of these eigenvalue problems, phase diagrams 
representing in the GS parameter space corresponding to the onset of spot instabilities are obtained for various simple 
spatial configurations of multi-spot patterns. In addition, it is shown that there is a wide parameter range where a spot 
instability can be triggered only as a result of the intrinsic slow motion of the collection of spots. The construction of 
the quasi-equilibrium multi-spot patterns and the numerical study of the spectrum of the eigenvalue problems relies 
on certain detailed properties of the reduced-wave Green's function. The hybrid asymptotic-numerical results for spot 
dynamics and spot instabilities are validated from full numerical results computed from the GS model for various spatial 
configurations of spots. 
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1 Introduction 

Spatially localized spot patterns have been observed in a wide variety of experimental settings including, the 
ferrocyanide-iodate-sulphite reaction (cf. |40j ). the chloride-dioxide-malonic acid reaction (cf. |14j l. and certain 
electronic gas-discharge systems (cf. [3], [4]). Furthermore, numerical simulations of certain simple reaction-diffusion 
systems, such as the two-component Gray-Scott model (cf. |48j . [42]) and a three-component gas-discharge model 
(cf. [52]), have shown the occurrence of complex spatio-temporal localized spot patterns exhibiting a wide range 
of different instabilities including, spot oscillation, spot annihilation, and spot self-replication behavior. A survey of 
experimental and theoretical studies, through reaction-diffusion (RD) modeling, of localized spot patterns in various 
physical or chemical contexts is given in [55] . These experimental and numerical studies have provided considerable 
impetus for developing a theoretical understanding of the dynamics, stability, and spot self-replication behavior asso- 
ciated with localized solutions to singularly perturbed RD systems. A brief survey of some new directions and open 
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problems for the theoretical study of localized patterns in various applications is given in [30] . More generally, a wide 
range of topics in the analysis of far-from-equilibrium patterns modeled by PDE systems are discussed in |45j . 

In this context, the goal of this paper is to provide a detailed case study of the dynamics, stability, and self- 
replication behavior, of localized multi-spot patterns for the Gray-Scott (GS) RD system in a two-dimensional 
domain. The GS model is formulated in dimensionless form as (cf. [44] ) 

= Aw - w + yluw^ , TUt = D /\u+ {I- u) - uv^ , x G 17 G ; dnU = dnV = G , x € dn. (1.1) 

Here £ > 0, Z? > 0, r > 1, and A are constants. The parameter A is called the feed-rate parameter. For various ranges 
of these parameters, (jl.ip has a rich solution structure consisting of oscillating spots, spot annihilation behavior, spot 
self- replication, and stripe and labyrinthian patterns. For the specific choice D = 2e^, the complexity and diversity 
of these patterns were first studied numerically in [48] in the unit square. Further numerical studies that clearly 
exhibit the distinguishing phenomena of spot self-replication in a two-dimensional domain include 42. and [50j . 

For the study of localized patterns in the GS model (jl.ip . there are two distinguished limits for the diffusion 
coefficient D in (jl.ip : the weak interaction regime with D = ©(e^), where the original numerical simulations of 
the GS model were performed (cf. [48]), and the semi-strong interaction regime D = 0(1), where many analytical 
studies have been focused. In the weak interaction regime, there is only an exponentially weak coupling between any 
two spots in the multi-spot pattern. This weak coupling arises from the exponential decay of a local spot profile. In 
contrast, in the semi-strong interaction regime, and for a certain range of A, the spots are more strongly coupled 
through the long-range effect of the slowly varying u component in (jl.ip . In this way, for D — 0(1) the dynamics of 
each individual spot in a multi-spot pattern is rather strongly influenced by the locations of the other spots in the 
pattern, as well as by the geometry of the confining domain. 

In a one-dimensional spatial domain, there has been much work over the past decade in analyzing the stability, 
dynamics, and self-replication of spike patterns for the GS model (|l.ip . For the weak interaction regime where 
D — O(e^) and A = 0(1), the mechanism for spike self-replication put forth in [46] (see also [54] ) was based on 
the occurrence of a nearly-coinciding hierarchical saddle-node global bifurcation structure for the global bifurcation 
branches of multi-spike solutions. This mechanism was also shown to occur for the related Gierer-Meinhardt (GM) 
system (cf. [45] ). In this one-dimensional context, it was shown in [24^ that typically only the spikes at the edges 
of a multi-spike pattern can undergo splitting. The possibility of spatial-temporal chaotic behavior of spike patterns 
due to repeated annihilation and self-replication events was explored in [47j from a global bifurcation viewpoint. The 
study of solution behavior in this weak interaction regime relies heavily on the use of numerically computed global 
bifurcation diagrams, since it appears to be analytically intractable to study the local problem near each spike. In 
contrast, for the semi-strong interaction regime D = 0(1) there are many analytical studies of spike behavior for 
(jl.ip for different ranges of the parameter A. For the range 0(e^/^) < A 0(1), oscillatory instabilities of the spike 
profile, characterized in terms of threshold values of the time-constant r in (jl.ip . have been analyzed in [17] . [18] . 
[43] . [34] . and [9]. Competition or annihilation instabilities of the spike profile, characterized by threshold values of 
the diffusivity D, have been analyzed in 34 for the range A = 0(e^/^). In addition, self-replication instabilities of 
spike patterns have been shown to occur only in the regime A = 0(1), and they have been well-studied in [50] , [51] , 
[18] . [44], [35] , and [20] . Weak translation, or drift, instabilities of spike patterns have been analyzed in [33^ and 
[35] . Finally, there have been several studies of the dynamical behavior of spike patterns for the one-dimensional GS 
model including, two-spike dynamics for the infinite-line problem (cf. |15] , |16) ) and in a bounded domain (cf. [53]), 
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and multi-spike patterns on a bounded domain (cf. |9]). Related studies on the stability and dynamics of spike 
solutions for the GM model in a one-dimensional domain are given in |29j . |28j . |59j . |21j . and [56j (see also the 
references therein). For the semi-strong regime D = one key feature of the GS model (|l.ip in one spatial 

dimension is that the parameter regime A = 0{e^^'^) where spike competition instabilities occur is well-separated in 
parameter space from the range A — 0(1) where spike self-replication occurs. As we discuss below, this feature with 
the one-dimensional GS model is in distinct contrast to the two-dimensional GS model (jl.ip where several distinct 
spot instability mechanisms occur in nearby parameter regimes for A. 

Although the stability properties of spike patterns for the one-dimensional case is rather well-understood, there 
are only a few studies of the stability of multi-spot patterns for singularly perturbed two-component RD systems in 
two dimensional domains. In particular, for the GS model (11.11) on the infinite plane 17 = M^, the existence and the 
stability, with respect to locally radially symmetric perturbations, of a one-spot solution to (II. ip was studied in [60) 
for the range A = 0{e{—\ney^'^) with either D = 0(1) or D = 0(j^^^), where = — 1/lne . This rigorous study 
was based on first deriving, and then analyzing, a certain nonlocal eigenvalue problem (NLEP). For the same range 
of A, in |62j the one-spot NLEP stability analysis of [60j was extended to treat the case of multi-spot patterns on 
a bounded domain. A further extension of this theory to study certain asymmetric multi-spot patterns was made 
in |61j . The A:-spot NLEP stability analysis of |62j for the regime A = 0(£(— Ine)^/^) was based on retaining only 
the leading-order term in powers of v = — 1/lne in the construction of the spot profile, and it pertains to the class 
of locally radially symmetric perturbations near each spot. This theory characterizes competition and oscillatory 
profile instabilities for the parameter range A = 0(e(-lne)i/2) with either D = 0(1) or D = 0{iy-^). In this 
leading-order theory, the stability thresholds depend only on the number of spots, and not on their spatial locations 
in the multi-spot pattern. Spot self-replication instabilities were not accounted for in these NLEP studies, as this 
instability is triggered by a locally non-radially symmetric perturbation near each spot for the nearby parameter 
regime A = 0{e{— Ine)) (cf. [43 ). A survey of NLEP stability theory as applied to other two-component singularly 
perturbed RD systems, such as the GM and Schnakenburg models, is given in [63]. 

With regards to the dynamics of spots, there are only a few analytical results characterizing spot dynamics for 
singularly perturbed RD systems in two-dimensional domains. These include, [31] and |58| for a one-spot 

solution of the GM model, [22], |23j and [25) for exponentially weakly interacting metastable spots in various 
contexts, |38] for the Schnakenburg model, and [52 for a three-component gas-discharge RD model. We are not 
aware of any previous study of the dynamics of spots for the GS model (jl.ll) in a two-dimensional domain. 

We emphasize that the previous NLEP stability studies for the GS model ()l.ip (cf. [60j . [61j . |62j ) are based on 
a leading-order theory in powers of v — — 1/ Ine for the parameter range A = 0{e{— Ine)^/^) with either D = 0(1) 
or D = 0(j^^^). Therefore, since v is not very small unless e is extremely small, it is desirable to obtain a stability 
theory for multi-spot solutions that accounts for all terms in powers of v. However, more importantly, since the 
scaling regime A — 0(— elne) where a spot- replication instability can occur for (|l.ip (cf. |43| ) is logarithmically 
close to the low feed-rate regime A — 0(e(— Ine)^/^) studied in [62 and |60j . where only competition or oscillatory 
profile instabilities can occur, it is highly desirable to develop an asymptotic theory that incorporates these two 
slightly different scaling regimes into a single parameter regime where all three types of spot instability can be 
studied simultaneously. The leading-order theory in [62] is not sufficiently accurate to study the three types of spot 
profile instability (competition, oscillatory, and self-replication) within a single parameter regime. 

For the simpler Schnakenburg RD model, where competition and oscillatory instabilities do not occur when D = 
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0(1), a hybrid asymptotic-numerical method was developed in [38] to study the dynamics and self-replication 
instabilities of a collection of spots for this specific RD model in an arbitrary two-dimensional domain. The theory, 
which accounts for all terms in powers oi v — —l/h\e, was illustrated explicitly for the square and the disk, and the 
results from this theory were favorably compared with full numerical computations of the RD system. 

The main goal of this paper is to extend the theoretical framework of [38] for the Schnakenburg model to study 
the dynamics and three types of instabilities associated with a collection of spots for the GS model in the 

semi-strong parameter regime D = 0(1) with A = ©(—sine). In our theory we account for all terms in powers of 
V = — 1/lne. In contrast to the Schnakenburg model of |38| . we emphasize that there are three distinct instability 
mechanisms for a collection of spots to the GS model in this scaling regime for A and D that must be considered. 

We now give an outline of the paper. In Sj2]and fj3]the method of matched asymptotic expansions is used to construct 
a quasi-equilibrium fc-spot pattern for (jl.ip that evolves slowly over a long 0{e~'^) time scale. For this pattern, the 
spatial profile for v concentrates at a set of points Xj S f2 for j = 1, . . . , fc that drift with an asymptotically small 
0{£^) speed. Within an 0{£) neighborhood of each spot centered at Xj, and at any instant in t, the local spot profiles 
for u and v are radially symmetric to within 0(e) terms and satisfy a coupled system of BVP, referred to as the core 
problem, on the (stretched) infinite plane. In the outer region, each spot at a given instant in time is represented 
as a Coulomb singularity for u of strength Sj. The collective coordinates characterizing the slow dynamics of this 
quasi-equilibrium /c-spot pattern are the locations xi , . . . , x^ of the spots and their corresponding source strengths 
Si, . . . , Sk, which measure the far-field logarithmic growth of the (inner) core solution for u near each spot. By 
asymptotically matching the inner and outer solutions for u, we derive a differential algebraic system (DAE) of 
ode's for the slow time evolution of these collective coordinates. At any instant in time, the quasi-equilibrium 
solution is characterized as in Principal Result 2.1, where the source strengths Sj for j — l,...,fc are shown to 
satisfy a coupled nonlinear algebraic system that depends on the instantaneous spot locations x^ for j = 1, . . . , fc, 
together with certain properties of the reduced- wave Green's function G'(x;Xj) and its regular part Rjj defined by 

AG- -^G= -(5(x-Xj), xer2; a„G = , x£dn, (1.2 a) 

G(x;x,-) ^ — In |x — xJ + .j + 0(1) , as x — > x,- . (1-2 &) 

27r 

In Principal Result 3.1 of SQlthe dynamical behavior of the collection of spots is characterized in terms of the source 
strengths and certain gradients of the reduced- wave Green's function. The overall DAE ODE system for Xj and Sj, 
for j = 1, . . . , fc, incorporates the interaction between the spots and the geometry of the domain, as mediated by 
the reduced-wave Green's function and its regular part, and it also accounts for all logarithmic correction terms in 
powers oi V ^ — 1 / In e in the asymptotic expansion of the solution. In this DAE ODE system there are two nonlinear 
functions of Sj , defined in terms of the solution to the core problem, that must be computed numerically. 

In ij4.1l we study spot self-replication instabilities by first deriving a local eigenvalue problem near the spot that 
characterizes any instability due to a non-radially symmetric local deformation of the spot profile. We emphasize that 
in this stability analysis the local eigenvalue problems near each spot are not coupled together, except in the sense 
that the source strengths Si, . . . , Sk must be determined from a globally coupled nonlinear algebraic system. The 
spectrum of the local two-component linear eigenvalue problem near the spot is studied numerically, and we show 
that there is a critical value S2 ~ 4.31 of the source strength Sj for which there is a peanut-splitting linear instability 
for any Sj > S2. These results for spot-splitting parallel those for the Schnakenburg model, as given in [38j . As a new 
result, we derive and then numerically study a certain time-dependent elliptic-parabolic core problem near the j^"'^ 
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spot. Our computations from this time-dependent core problem strongly suggest that the localized peanut-splitting 
linear instability of the quasi-equilibrium spot profile is in fact subcritical, and robustly triggers a nonlinear spot self- 
replication event for the spot whenever Sj > S2. In summary, our hybrid asymptotic-numerical theory predicts 
that if S'j > S2 ~ 4.31 for some J € {1, . . . ,k} then the spot will undergo a nonlinear spot self-replication event. 
Alternatively, the spots are all stable to self-replication whenever 5j < S2 for j = 1 . . . , fc. 

In ii4.2l we use the method of matched asymptotic expansions to formulate a novel global eigenvalue problem 
associated with competition or oscillatory instabilities in the spot amplitudes for a fc-spot quasi-equilibrium solution 
to (|l.ip for the parameter range A = 0{—elne) with D = 0(1). This global eigenvalue problem, as formulated in 
Principal Result 4.1, is associated with a locally radially symmetric perturbation near each spot, and it accounts 
for all terms in powers of i^. It differs from the eigenvalue problem characterizing spot self-replication in that now 
the local eigenfunction for the perturbation of the u component in p.ip has a logarithmic growth away from the 
center of each spot. This logarithmic growth leads to a global eigenvalue problem that effectively couples together 
the local problems near each spot. A key component in the formulation of this global eigenvalue problem is a certain 
eigenvalue-dependent Green's matrix, with entries determined in terms of properties the Green's function G'a(x;Xj) 
satisfying AG — D^^{1 + t\)G = —S{x — Xj) for x G fi, with 9„G = for x e dfl. This globally coupled eigenvalue 
problem can be viewed, essentially, as an extended NLEP theory that accounts for all terms in powers of ly. In 
^5]-^ we show that it determines thresholds for competition and oscillatory instabilities very accurately. However, in 
contrast to the leading-order-in-j/ NLEP stability studies (cf. |60j . [62]) for A~0 {^e[— Ine]^/^), our globally coupled 
eigenvalue problem is not readily amenable to rigorous analysis. In Appendix B we briefly review the NLEP theory 
of |60j and |62j . and we show how our globally coupled eigenvalue problem can be reduced to leading order in v to 
the NLEP problems of gS] and |62] when A^O Ine]!/^) and D = 0{v-^). 

In our stability analysis of we linearize the GS model (jl.ip around a quasi-equilibrium solution where the spots 
are assumed to be at fixed locations Xi, . . . ,Xfe, independent of time. However, since the spots locations undergo 
a slow drift with speed ©(e^), the source strengths Sj for j — 1, . . . ,fc also vary slowly in time on a time-scale 
t — ©(e^^). As a result of this slow drift, there can be triggered, or dynamically induced, instabilities of a quasi- 
equilibrium spot pattern that is initially stable at time t = 0. To illustrate this, suppose that the pattern is initially 
stable to spot self-replication a± t — Q in the sense that Sj < S2 at t = for j 1, . . . , fc. Then, it is possible, 
that as the spot drifts toward its equilibrium location in the domain, that Sj > E2 after a sufficiently long 
time of order t = 0{s^'^). This will trigger a nonlinear spot self-replication event for the J^^ spot. In a similar way, 
we show that dynamically-triggered oscillatory and competition instabilities can also occur for a multi-spot pattern. 
This dynamical bifurcation phenomena is similar to that for other ODE and PDE slow passage problems (cf. [5], 
[39j ) that have triggered instabilities generated by a slowly varying external bifurcation, or control, parameter. The 
key difference here, is that the dynamically-triggered instabilities for the GS model p.ip occur as a result of the 
intrinsic motion of the collection of spots, and is not due to the tuning of an external control parameter. 

In our numerical computations of competition and oscillatory instability thresholds from our globally coupled 
eigenvalue problem of S21 we will for simplicity only consider fc-spot quasi-equilibrium spot configurations xi , . . . , 
for which a certain Green's matrix is circulant symmetric. For instance, this circulant matrix structure occurs when 
k spots are equally spaced on a circular ring that is concentric within a circular disk, and it also occurs for other 
spot patterns with sufficient spatial symmetry in other domains. Examples of such patterns are given in [J5] and |6] 
below. Under this condition, we show in t ^2.1l that the source strengths Sj for j — 1, . . . , fc have a common value. In 
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addition, by calculating the spectrum of the circulant symmetric Green's matrix, we show in Principal Result 4.3 of 
§4.3l that the globally coupled eigenvalue problem simplifies to k separate transcendental equations for the eigenvalue 
parameter. In Appendix C we outline the numerical methods that we use to compute the instability thresholds from 
the globally coupled eigenvalue problem under the circulant Green's matrix assumption. 

Spot patterns that give rise to this special circulant matrix structure are the direct counterpart of equally- spaced k- 
spike patterns with spikes of a common amplitude, treated in almost all of the previous NLEP stability studies of the 
GS and related RD models on a one-dimensional domain (cf. [l7j, [18], [34], [15], [16], [29], [59], [21], [56]). In one 
spatial dimension, the only NLEP stability studies of an arbitrarily- spaced slowly evolving fc-spike quasi-equilibrium 
solution are the asymptotic-numerical study of dynamic competition instabilities for the Gierer-Meinhardt model 
with r = in 28 , and the study of oscillatory instabilities in [9] for the one-dimensional GS model for the range 
0(61/2)^^^0(1). For this range of A it was shown in [9] that the k separate NLEP problems can be reduced, via 
a scaling law, to only one single NLEP problem. To date, there has been no NLEP stability study of a slowly evolving 
arbitrarily- spaced fc-spike quasi-equilibrium spike pattern in a one-dimensional domain that takes into account both 
competition and oscillatory instabilities. As a result, in our two-dimensional setting, it is a natural first step to study 
the global eigenvalue problem, which governs competition and oscillatory instabilities, under the circulant matrix 
condition, which allows for spots of a common source strength. 

In ^the asymptotic theory of Sj2]- 21 is illustrated for the case of both one and two-spot quasi-equilibrium solutions 
to the GS model (II. ip on the infinite plane. Phase diagrams characterizing the GS parameter ranges for the different 
types of instabilities are derived for these simple spot patterns. In particular, for two spots that are sufficiently 
far apart, we show that spot self-replication instabilities will occur when A exceeds some threshold. In contrast, a 
competition instability will occur if the two spots are too closely spaced. For a very large, but finite, domain the full 
numerical simulations in H5.3I are used to validate the stability results from the asymptotic theory. 

In Sj6] the asymptotic theory of 21 is implemented and compared with full numerical results computed from 
p.ip for various special multi-spot patterns on the unit disk and square for which a certain Green's matrix has 
a circulant matrix structure. For these domains, the explicit formulae for the reduced-wave Green's function and 
its regular part, as derived in Appendix A, are used to numerically implement the asymptotic theory. The overall 
hybrid asymptotic-numerical approach provides phase diagrams in parameter space characterizing both the stability 
thresholds and the possibility of dynamically-triggered instabilities. One key theoretical advantage of considering the 
unit disk is that the reduced-wave Green's function can be well-approximated for 13 ^ 1 by the Neumann Green's 
function, which has a simple explicit formula in the unit disk. By using this explicit formula, the hybrid asymptotic- 
numerical framework of oilcan be studied, to a large extent, analytically for the case of k equally-spaced spots 
on a ring that is concentric within the unit disk. 

In ^ we compare our theoretical predictions for spot dynamics and spot self-replication instabilities with full 
numerical results computed from (|l.ip for a few simple "asymmetric" spot patterns for which the associated Green's 
matrix is not circulant. We show that the dynamics in Principal Result 3.1 accurately determines spot dynamics before 
a self-replication event, and with a re-calibration of the initial spot locations, it accurately predict spot dynamics 
after a spot-splitting event. We emphasize that since the local eigenvalue problems near each spot are decoupled for 
the case of locally non-radially symmetric perturbations, the onset of spot self-replication behavior only depends on 
the source strength of an individual spot. Therefore, given any initial spatial configuration of spots at t = 0, we need 
only solve the nonlinear algebraic system for Sj,i = l,--- , fc at t = to predict that the spot undergoes splitting 
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starting at t = when Sj > S2 ~ 4.31. A special asymmetric pattern that we consider in some detail in ^J7] is a 
/c-spot pattern consisting of A: — 1 equally-spaced spots on a ring concentric within the unit disk, with an additional 
spot at the center of the unit disk. For D ^ 1, we use the simple explicit formula for the Neumann Green's function 
to explicitly predict the occurrence of dynamically-triggered spot self-replication instabilities for this special pattern. 

Although the hybrid asymptotic-numerical framework developed herein to study the stability and dynamics of 
multi-spot quasi-equilibrium patterns for the GS model (jl.ip is related to that initiated for the Schnakenburg model 
in [38] ■ there are some key differences in the analysis and in the results obtained. The primary difference between 
these two models is that the GS model (|l.ip admits three types of instability mechanisms, whereas only spot self- 
replication instabilities can occur for the Schnakenburg model of |38j when D = 0(1). In addition, in contrast to 
our study in ^JSlof one- and two-spot patterns to the GS model on the infinite plane, the Schnakenburg model of [38] 
is ill-posed in R^. Finally, our results for the GS model (|1.1|) show that there is a wide parameter range and many 
simple spot configurations for which we can theoretically predict the occurrence of dynamically-triggered instabilities 
due to either competition, oscillation, or splitting. These dynamically-triggered bifurcation events can occur even 
within the very simple context of a multi-spot pattern with a common source strength. For these special patterns, 
such instabilities cannot occur for the Schnakenburg model. 

Finally, although this paper focuses only on the study of spot patterns, we remark that the GS model ()1.1|) supports 
patterns of increasing complexity as the feed-rate parameter A increases. In particular, in the range 0{e^^^) <C A <C 
0(1), the GS model (II. ip with D = 0(1) on a two-dimensional domain does not admit spots, but instead allows 
for solutions for which v concentrates on a higher dimensional set such as on a one-dimensional stripe or a one- 
dimensional ring inside a two-dimensional domain. A stability analysis of a planar stripe inside a square domain or 
a concentric ring inside a disk was given in j41] and in [37j . 



2 if-Spot Quasi-Equilibrium Solutions 

We first construct a fc-spot quasi-equilibrium solution to (jl.ip by using the method of matched asymptotic expansions. 
We denote the center of the spot by Xj — {xj, i/j) G 17 for j — 1, . . . , fc. We assume that the spots are well-separated 
in the sense that |xi — Xj| ~ 0(1) for i 7^ j, and dist(xj, dfl) — 0(1) for j = 1, . . . , k. In an 0(e) neighborhood near 
the j^^ spot, we get v = 0{e^^) and u — 0(e). Thus, we introduce the local variables Uj, Vj, and y, defined by 



^ Uj, v^—V,, y = e-i(x-x,). (2.1) 



AVD - e 
In terms of these local variables, ()l.ip transforms to 

^yVj - V, + U,Vf - , Ayt/, - U,Vf + - f!c/^. = , y G (2.2) 



We look for a radially symmetric solution to (12.21) of the form IJj = Vji^p) and Vj = ^^(p), where p= \y\. Then, to 
leading order in e, \] j and Vj are the solutions to the radially symmetric problem 

C/j' + ^C/j - t/jF/ = 0, V'j + ^Vl -Vj+UjVf = ^ , 0<p<oo, (2.3 a) 

F/(0) = , [/j (0) = ; V, (p) ^ , L/j (p) - In p + x{Sj) + o(l) , as p ^ 00 . (2.3 h) 

This leading-order coupled inner problem is referred to as the core problem, and is the same as that derived in [38] 
for the Schnakenburg model. We refer to Sj as the source strength of the spot. From the divergence theorem, it 
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follows from the Uj equation in (|2.3p that Sj — UjV^pdp > 0. In the far-field behavior (|2.3 for Uj, the constant 
X is a nonlinear function of the source strength Sj , which must be computed numerically from the solution to ()2.3|) . 

The solution to (|2.3p is calculated numerically for a range of values of Sj > by using the BVP solver COLSYS 
(cf. ^). In Fig. m we plot xi^j), ^(0) versus Sj, and Vj{p) for a few different values of Sj. For Sj > Sy ^ 4.78, 
the profile Vj{p) has a volcano shape, whereby the maximum of Vj occurs at some p > 0. These computations give 
numerical evidence to support the conjecture that there is a unique solution to (|2.3p for each Sj > 0. 






(a) X vs. Sj 



(b) 1/ (0) vs. S, 



(c) Vj vs. p 



Figure 1. Numerical results for the core problem (j2.3|) (^a^ The function x 'vs. Sj; (b) Vj{0) vs. Sj; (c) The spot profile 
Vj{p) for Sj = 0.94, 1.45, 2.79 (solid curves where Vj(0) is the maximum ofVj{p)), and the volcano profile Vj{p) for 
Sj = 4.79, 5.73, 6.29 (the dotted curves correspond to cases where the maximum ofVj{p) occurs for p > 0). 

Since v is localized near each Xj for j ~ 1, . . . , fc, and is exponentially small in the outer region away from the spot 
centers, the effect of the nonlinear term uv'^ in the outer region can be calculated in the sense of distributions as 
uv'^ ^ Y^j=\ ( /m2 (^s)"^ ^o^f ^y) '^(^ ~ ^j) ^ 2Tre^/DA~^ Sj=i ^^j)- Therefore, in the quasi-steady 

limit, the outer problem for u from (jl.ip is 



DAu + (1 - u) = ^""^^^ ^ Y Sj (5(x ~ Xj) , xer!; 



A 



dnU = , X e 9r2 , 



In |x — Xj I — Sj Ine + x{Sj)j , as x — , j = 1, . . . , fc . 



(2.4 a) 



(2.4 &) 



The singularity condition (j2.4 6P for u as x — > Xj was derived by matching the outer solution for u to the far-field 
behavior (j2.3 b\i of the core solution, and by recalling u = ellj / (A^/D) from (|2.ip . The problem (|2.4I) suggests that 
we introduce new variables A = 0(1) and v ^ 1 defined by 

i/ = -l/lne, y^ = zyAV!D/£ = A/D/[-eln£]. (2.5) 

In terms of these new variables, (12.41) transforms to 

fc 

(2.6 a) 



Am - 



(1 



9„M = , X g 9ri , 



1 



[S'jj^ln |x - Xjl + Sj + i^x{Sj)] 



as X x," 



j = l,...,k. 



(2.6 6) 



We emphasize that the singularity behavior in (|2.6 6[) specifies both the strength of the logarithmic singularity for u 
and the regular, or non-singular, part of this behavior. This pre-specification of the regular part of this singularity 
behavior at each Xj will yield a nonlinear algebraic system for the source strengths Si, . . . , Sk. 
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The solution to (12. 6|) is represented as u = 1 — X^iLi ^ni^A^^ SiG{yi; x^), where G(x; x^) is the reduced- wave Green's 
function defined by ()1.2p . By expanding u as x — > Xj, and then equating the resulting expression with the required 
singularity behavior in (j2.6 &|) . we obtain the following nonlinear algebraic system for Si, . . . , 5*^: 

k 

^=5'j(l + 27ri.i?j-j) + i^x(S'j) + 2vri/^S',G(xj;x,), j^l,...,k. (2.7) 

i=l 

Given the GS parameters A, e and D, we first calculate A and v from ()2.5p . and then solve ()2.7p numerically for the 
source strengths Si, . . . , Sk- With 5j known, the quasi-equilibrium solution in each inner region is determined from 
()2.3p . As a remark, since A = 0{l), then A — ©(— elne) from (12. 5p . Therefore, the error made in approximating 
(12. 2p by (|2.3p in the inner region is of the order ©(— Ine). We summarize our result as follows: 

Principal Result 2.1: For e ^ assume that A — 0{—e\ne), and define v and A hy v = — 1/lne and A = 
vA^/D/e, where A = 0{\). Then, the solution v and the outer solution for u, corresponding to a k-spot quasi- 
equilibrium solution of the GS model (jl.ip . are given asymptotically by 

k I — h 

«(x) ^ 1 - ^ ^ 5,G(x; X,) , «(x) - ^ E (^"'l'^ - I) • (2.8 a) 

Moreover, the inner solution for u, defined in an 0{£) neighborhood of the j^^ spot, is 

u(x)^^C/,(£-i|x-x,|) . (2.8 6) 

Here xi,...,Xfe is the spatial configuration of the centers of the spots, and G(x;Xj) is the reduced-wave Green's 
function satisfying ()1.2p . In (j2.8p . each spot profile Vj{p) and Uj{p) for j = 1, . . . , fc satisfies the coupled EVP system 
(|2.3p . where the source strength Sj in (j2.3 &P is to be calculated from the nonlinear algebraic system (|2.7p . 

We emphasize that the nonlinear algebraic system (12. 7p determines the source strengths Sj for j = 1 , . . . , fc to 
within an error smaller than any power of z/ = —1/ Ine. As such, our construction of the quasi-equilibrium pattern 
is accurate to all orders in v. Similar techniques for summing logarithmic expansions in the context of linear elliptic 
PDE's or eigenvalue problems in two-dimensional perforated domains containing small holes have been developed 
in a variety of contexts (cf. [57j . |32] . |13) . [49j). Our construction here of a quasi-equilibrium solution for the GS 
model (jl.ip extends this previous methodology for treating logarithmic expansions to an RD system for which the 
local problem is nonlinear with a logarithmic far-field behavior. The nonlinear algebraic system (|2.7p for the source 
strengths is the mechanism through which the spots interact and sense the presence of the domain fl. This global 
coupling mechanism, which is not of nearest-neighbor type as in the case of the exponentially weak spot interactions 
studied in 1231 and [25|, is rather significant since v = — 1/lne is not very small unless e is extremely small. 

The quasi-equilibrium solution in Principal Result 2.1 exists only when the spatial configuration of spots and the 
GS parameters are such that the nonlinear algebraic system (|2.7p for Si, . . . , Sk has a solution. Determining precise 
conditions for the solvability of this system is a difficult issue. Therefore, in ^JSHZ] we will primarily consider spot 
patterns where the source strengths have a common value, for which (j2.7p reduces to a scalar nonlinear equation. 

Next, we derive analytical approximations for the solution to (|2.7p by first re-writing (j2.7p in matrix form. To do 
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so, we define the Green's matrix the vector of source strengths s, the vector x(s), and the identity vector e by 
/ Gi^2 • • • Gi^k \ 

^2,1 ■• ■• : 

Sk I \ 1 



Q = 



Si 



1 



/ xiSi) 



e = 



X(s) 



(2.9) 



\ xiSk) 



■ ■• ■• Gk-i,k 

V Gk,l ■ ■ ■ Gk,k-1 Rk,k / 

Here Gij = G{xi;xj), and Gij — Gj^i by reciprocity, so that Q is a. symmetric matrix. Then, p.7p becomes 



Ae = s + 2Tnygs + iyx{s). (2.10) 

We will consider (|2.10p for two ranges oi D; D ^ £>(!) and D = 0{i'^^). For D = 0(1), we can obtain a two-term 
approximation for the source strengths in terms of ^ 1 by expanding s = So + i^Si + • • • . This readily yields 

s = Ae-i' {2-kA a e + x{A)e\ + 0{v^) . (2.11) 

Therefore, for v ^ \, the leading-order approximation for s is the same for all of the spots. However, the 0(y) 
correction term depends on the spot locations and the domain geometry. 

Next, we consider (|2.10l) for the distinguished limit where D = Do/iy ^ 1 with Do = 0(1). Since the reduced-wave 
Green's function G, satisfying ()1.2|) . depends on D we first must approximate it for D large. Assuming that f2 is a 
bounded domain, we expand G and its regular part R for _D 3> 1 as 

G - DG-i + Go + -^Gi + • • • , i? - Z3i?_i + i?o + ^Ri + ■■■ ■ 

Substituting this expansion into (jl.2[) . and collecting powers of D, we obtain that G_i is constant and that 

AGo = G_i -(5(x-xj), xeQ: c»„Go = , x e dn . 

The divergence theorem then shows that G_i = where is the area of fl. In addition, the divergence theorem 

imposed on the Gi problem enforces that Gq dx = 0, which makes Go unique. Next, from the singularity condition 



(|1.2 6p for G we obtain for x — >■ Xj that DG-i + Go(x;xj) 



-(27r)-iln|x-Xj| + DR^i +i?o(x;Xj) + 



Since G-i = ^, we conclude that = G-i = l/|i^|. In this way, we obtain the following two-term expansion 

G(x;x,)^G(x;x,) = ^ + GW(x;x,) + ... , for I? » 1 . (2.12) 

Here GW(x;Xj) is the Neumann Green's function with regular part -R^-^\ determined from the unique solution to 



R 



'3,J 



- \n\ 



R 



(N) 
3,j 



AG(^) = -^-(5(x-Xj), xeQ: dnG^^'^:^0, x e 9f] ; / G^^^ dx = , 
l^^l Jn 

^(-)-o(l). 



GW(x;x,)^-^ln|x-x, 



as X 



(2.13 a) 
(2.13 6) 



For the disk and the square, in Appendix A we analytically calculate both G(x;xj) and its regular part Rjj, as 
well as the Neumann Green's function G^^'' (x; Xj ) and its regular part -R^^'' • For the case of a one-spot solution 
centered at the midpoint xi of the unit square [0, 1] x [0, 1], we use some of the explicit formulae from Appendix A 
to compare the two-term approximation given in (j2.12p . with the reduced- wave regular part as computed 
from (jA.lip . These results are shown in Fig. [2] From (jA.13p . the two-term approximation p.l2p for large D is 



^1, 



D 



12 27r 



ln(27r) , 



-27r 



(2.14) 
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A similar comparison is made in Fig. [3]for the case of a single spot located at the center xi — (0, 0) of the unit disk. 
For this radially symmetric case Ri,i can be found explicitly, and its two-term approximation for large D is obtained 
from (|2.12p and (|A.9 &[) . In this way, we get 



ilni:i + ln2-7e 



rii.i — . 

IT OTT 



(2.15) 



Here 7e ~ 0.5772 is Euler's constant, while Ii{r) and Ki{r) are the modified Bessel functions of order one. From 
Fig. [Hand Fig. [3] we observe that for both domains the two-term approximation (|2.12l) involving the regular part of 
the Neumann Green's provides a decent approximation of Ri^i even for only moderately large values of D. 




(a) and vs. D 




(b) Ri.i and vs. D 



Figure 2. Consider a single spot centered at the midpoint xi = (0.5, 0.5) of the unit square. We plot Ri i vs. D (solid 
curve) and its two-term large D approximation given in p.l4p (dotted curve), (a) D e [0.1,6]; (b) D e [0.1,0.5]. 



1.5 

1 

0.5 


-0.5 



DC 



12 3 4 

D 

(a) and Ri^i vs. D 




(b) Ri.i and vs. D 



Figure 3. Consider a single spot located at the center Xi = (0, 0) of the unit disk. We plot Ri i vs. D given in 
i2.15\) (dotted curve) and its two-term large D approximation given in 112.15]) (solid curve), (a) D G [0.1,6]; (b) 
D e [0.1,0.5]. 

Next, we use the large D asymptotics to find an approximate solution to the nonlinear algebraic system l|2.10p in 
the hmit D = Dq/i/ > 1, where Dq = 0(1) and < 1. Upon substituting (|2.12|) into (|2.10|) we obtain 



= s 



^ee^s + 27r;.gWs + ^.x(s), 



where Q^^'^ is the Green's matrix associated with the Neumann Green's function, i.e. ^I^"* = G'^^^(xi;xj) for i ^ j, 
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and ^j^'' = ^jT'' ^■^ expanding s as s = Sq + i^si + • • • for <C 1, we then obtain that Sg and Si satisfy 

+ ^ee^) so ^ , {l + ^ee^) = -2^g^^)s, ~ xi-o) , (2.16) 

where / is the k x k identity matrix. Since e-^e — k, the leading-order approximation Sq shows that the source 
strengths have an asymptotically common value Sc given by 

The next order approximation Si from (j2.16l) yields 

si = - (/ + ^^ee^y^ ((Z^S^G^'''^ + x(5c))e . (2.18) 

Since the matrix / + fiee^ is a rank-one perturbation of the identity, its inverse is readily calculated from the 
Shermann- Woodbury-Morrison formula as (/ -I- fiee^)^^ — I ~ ^ee'^/(l -I- ^k), which determines Si from (|2.18p . In 
this way, for D — Dq/v and ^ 1 we obtain the two-term expansion for s given by 



1 + fik 1 + fik\ 1 + fi 

Here Sc and fi are defined in (j2.17p . while the scalar function i^(xi, . . . , x^) is defined by 



Fta ^.)=e-^g<">e = f2EeS'- (2-20) 

i=l J = l 

In contrast to the leading-order approximation in (12. lip when D = 0(1), the leading-order approximation Sc in 
()2.19p depends on the number of spots and the area of the domain, with Sc increasing as the area \Q\ increases. 

We now illustrate our asymptotic theory for the construction of quasi-equilibria for the case of a one-spot solution 
centered at the midpoint of cither the unit square or disk. Many additional examples of the theory are given in ^JSHH] 

For our first example, we consider a one-spot solution with a spot located at the center xi = of the unit disk with 
e = 0.02. We fix = 1, so that from (|2.15p the regular part of the Green's function for a spot at xi is « 0.1890. 
Then, (12. 7p reduces to the following scalar nonlinear algebraic equation for the source strength Si in terms of A: 

A = Si{l + 2TT,yRi^i) + i^x{Si). (2.21) 



In Fig. 4(a) we plot A versus Si, showing the existence of a fold point at Af ~ 2.55 corresponding to Sf ~ 1.01. 
Thus, A > Af is required for the existence of a one-spot quasi-equilibrium solution located at the center of the unit 
disk. In Fig. |4(b)| the asymptotic result for u(xi ) = vUi (0) /A vs. ^ is shown by the dotted curve, with a fold point at 
u/(xi) w 0.50. The fold point is marked by a circle in both figures, and the critical value Ay ~ 5.55, u^(xi) w 0.083 
for a volcano-type solution corresponding to Sv ~ 4.78 is marked by a square. For A> Af, u(xi) has two solution 
branches. The upper branch corresponds to Si < Sf, while the lower branch is for the range > Sf. 

To validate the asymptotic result for solution multiplicity, we solve the steady-state GS model (|l.ip in the unit 
disk by using the Matlab EVP solver BVP^C. By varying u(0), we then compute the corresponding value of A. The 
resulting full numerical result for w(0) vs. A is shown by the heavy solid curve in Fig 4(b) which essentially overlaps 
the asymptotic result. This shows that when e = 0.02, the asymptotic result for the bifurcation diagram, based on 
retaining all terms in powers of v, agrees very closely with the full numerical result. 

For our second example, we consider a one-spot solution centered at the midpoint xi = (0.5, 0.5) of the unit square 
= [0, 1] X [0, 1]. We fix e = 0.02 and D = 1. Then, by using (jA.llP for the reduced-wave Green's function, as given 
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012345678 ' ' ' ' ' 

e 3 4 5 6 7 

A 

(a) A vs. 5*1 (b) u(xi) vs. A 



Figure 4. Let = {x||x| < 1}, e = 0.02, D = 1. and xi = 0. (a) A vs. Si; the square marks the volcano threshold 
Sv ~ 4.78, and the circle marks the fold point A/ ~ 2.55 at which Sf « 1.01. (b) m(xi) vs. A; the square marks 
Sy « 4.78, and the circle marks the fold point Sf ~ 1.01. The upper branch is for Si < Sf, and the lower branch is 
for Si> Sf. 




A 

(a) A vs. Si (b) u{xi) vs. A 



Figure 5. Let fl = [0,1] x [0,1], e = 0.02, D = 1, and Xi = (0.5,0.5). (a) A vs. Si; the square marks the volcano 
threshold Sy ~ 4.78, and the circle marks the fold point Af ~ 3.3756 at which Sf « 0.7499. (b) m(xi) vs. A; the 
square marks Sy ~ 4.78, and the circle marks the fold point Sf sa 0.7499. The upper branch is for Si < Sf, and the 
lower branch is for Si > Sf. 

in Appendix A, we obtain that Ri^i ~ 0.7876. We remark that even if we use the two-term large D asymptotics 
(j2.14p for the unit square, then we get the rather good estimate Ri^i w 0.7914. In Fig. [5] we use (|2.2ip to plot A 
versus Si and u(xi) ~ h'Ui{0)/A versus A. For A > Af ~ 3.376, w(xi) versus A has two solution branches, with 
the upper branch corresponds to Si < S f ~ 0.750 and the lower branch for 5*1 > Sf. In contrast to the previous 
example for the unit disk, where the full GS model (|1.1|) is reduced to a coupled set of BVP for ODE's, we cannot 
readily verify this asymptotic result for solution multiplicity in a square from a full numerical solution of (|l.ip . 

2.1 Symmetric Spot Patterns and a Circulant Matrix 

A special case for multi-spot patterns, which features prominently in [jJHT] below, is when the spatial configuration 
Xi , . . . , Xfc of spots within D, is sufficiently symmetric so that the Green's matrix Q is a circulant matrix. 

When Q is circulant, then it has the eigenpair Qe = 9e, where 9 — ^""'^ (^)i j- this special case. 
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(|2.10p has a solution for which the spots have a common source strength Sj = Sc for j = 1, . . . , fc, where Sc is the 
solution to the single nonlinear algebraic equation 

A = Sc + 2TnyeSc + vx{Sc) ■ (2.22) 

For 1/ <C 1, a two-term approximation for this common source strength Sc is 

Sc = A-i'[2tt9A + x{A)]+0{i'^). (2.23) 

For the distinguished limit D = 0{i'^^), and with Q a circulant matrix, we calculate from (|2.19l) that a two-term 
asymptotic approximation for the source strengths is given in terms of i^(xi, . . . , x^), as defined in (j2.20p . by 

Here 6''-^-' = F/k is the eigenvalue for the eigenvector e of the circulant Neumann Green's matrix Q^-^\ 

Owing to the fact that the nonlinear algebraic system (|2.10p can be reduced to the scalar nonlinear problem (|2.22p 
when the Green's matrix is circulant, in the majority of our numerical experiments for multi-spot patterns in SjS-fjT] 
below we will consider fc-spot quasi-equilibrium patterns that lead to this special matrix structure. 



3 The Slow Dynamics of a Collection of Spots 

In this section we derive the slow dynamics for the spot locations corresponding to a fc-spot quasi-equilibrium solution 
of the GS model At each fixed time t, the spatial profile of the spot pattern is characterized as in Principal 

Result 2.1. In the inner region near the spot, we introduce y — e^^ [x — Xj(^)], where ^ = e^t is the slow time 
variable, and we expand the inner solution as 



w=-4=(C/o,(p)+eC/i,(y) + ...) , i;-^(^o,(p)+el4,(y) + ...) . (3.1) 

Here the subscript in C/ojj^j denotes the order of the expansion, while j denotes the inner region. In the 
analysis below we omit the subscript j if there is no confusion in the notation. The leading-order terms U^j and Vqj 
are solutions of the core problem (|2.3p . Define Wj = (V^ij, Uij)'^ , where T denotes transpose. At next order, we get 
from p.ip and (jl.ip that Wj satisfies 

In p.2p . • denotes dot product, and eg = {cos 9, sin^)^, where 6 is the polar angle for the vector (x — Xj). 

To determine the dynamics of the spots we must calculate the gradient terms in the local expansion of the outer 
solution for u, as given in (j2.8p . in the limit x — > x^. To do so, we must calculate a further term in the local behavior 
as x — > Xj of the reduced- wave Green's function satisfying ()1.2p . In terms of the inner variable y, we get 

1 

2^ 



G(x; x,) ^ - — In |y| + R,,j + eVi?(xj; x^) • y - 



where we have defined Vi?(xj; Xj) = Vxi?(x; Xj 



By comparing the higher-order terms in the matching condi- 



tion between the inner expansion p.ip for u and the outer expansion for u, we obtain the required far-field behavior 
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as |y| — > oo of the inner solution Uij. In this way, we obtain that wj satisfies (|3.2p with far-field behavior 



_i..y J ' as y ^ oo , = i^y'^ ) =2tt\ 5, Vi?(x,; x,) + ^ 5,VG(x,; x,) | . (3.3) 



i = l 



To determine the dynamics of the spot we must formulate the solvability condition for (|3.2I) and (I3.3p . We 
define Pj{p) — {4'j{p)T'4'j{p))'^ to be the radially symmetric solution of the adjoint problem 

ApP* + MjP* = , < p < oo , (3.4) 

subject to the far-field condition that P* (0, l/p)^ as p — > oo, where Ap = dpp-\- p^^dp — p^^ . We look for solutions 
Pj and Pj to the homogeneous adjoint problem AyP* -\- MjP* = for y e in the form P* — Pj = P* cos and 
P* = Pj = P* sm0, where Pj{p) is the radially symmetric solution of (I3.4p . 

In terms of the adjoint solution Pj, the solvability condition for (I3.2p . subject to p.3p . is that 



lim / Pj ■ gj dy = lim 



Pj ■ 9pWj - Wj • dpP^ 



dy. (3.5) 

p=(T 



Here is a ball of radius a, i.e |y| = a. Upon using the far-field condition p.3p . and writing Xj = {xji,Xj2)'^ in 
component form, we reduce (j3.5p to 

^'27r POO p2iT 



/ VqCOS^ 9 p dp de -x'j^ / <^*l/o'cos6lsine'pdpd0 = lim / ( ij-ee)(jde. (3.6) 

Jo Jo Jo Jo '^^°°Jo V f / 

Therefore, since J^^ cos6sm6dd = 0, we obtain dxji/dS^ — 2/ji/ ^ (i)*V^pdpj . Similarly, the solvability condition 

for p.2p . subject to p.3p . with respect to the homogeneous adjoint solution P-, yields dxj2/d^ = 2/^2/ ^ /q°° (j)*V^pdp^ . 
Upon recalling the definition of {fji, fj2)^ in ()3.3p . we can summarize our result for the slow dynamics as follows: 

Principal Result 3.1: Consider the GS model with £ <C 1, ^ = 0(— elne), and t ^ Oie^"^). Then, provided 

that each spot is stable to any profile instability, the slow dynamics of a collection Xi , . . . , Xfc of spots satisfies the 
differential- algebraic (DAE) system 

^ ^ ~2^e^l{S,) 5, Vi?(x,; X,) -f ^ 5,VG(x,; x,) ] , j = 1, . . . , fc ; ^{S,) ^ . (3.7) 



S,<l>*Vipdp 



In (|3.7p . the source strengths Sj, for j — 1, . . . are determined in terms of the instantaneous spot locations and 
the parameters A and v of (|2.5p by the nonlinear algebraic system (|2.7p . In the definition ofj(Sj), Vq satisfies the 
core problem (|2.3p . while (j)* is the first component of the solution to the radially symmetric adjoint problem p.4p . 
Finally, the equilibrium spot locations Xje and spot strengths Sje, for j — 1, . . . , fc, satisfy 

k 

5'jeVi?(Xje;Xje) = ^' J^^^---^^, (3.8) 

i=l 

subject to the nonlinear algebraic system (12.71) . which relates the source strengths to the spot locations. 

The ODE system p.7p coupled to the nonlinear algebraic system (12. 7p constitutes a DAE system for the time- 
dependent spot locations Xj and source strengths Sj for j = 1, . . . , k. These collective coordinates evolve slowly over 
a long time-scale of order t — ©(e^^), and characterizes the slow evolution of the quasi-equilibrium pattern. From 
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a numerical computation of the function 7(>5'j) in (|3.7p was previously computed numerically in Fig. 3 of |38j , 
where it was shown that 7(5*^) > when Sj > 0. This plot is reproduced below in Fig. |10(b)[ In Sj7]we will compare 
the dynamics p.7p with corresponding full numerical results for different spot patterns in the unit square. 

We emphasize that the DAE system in Principal Result 3.1 for the slow spot evolution is only valid if each spot 
is stable to any spot profile instability that occurs on a fast 0(1) time-scale. One such spot profile instability is the 
peanut-splitting instability, studied below in ^14. 11 that is triggered whenever Sj > S2 ~ 4.31 for some J € {1, . . . , k}. 
The other profile instabilities are locally radially symmetric instabilities and, roughly speaking, consist of a temporal 
oscillation of the spot amplitude if t is sufficiently large, or a spot over-crowding competition instability, which is 
triggered when cither the spots are too closely spaced or, equivalently, when D is too large. A new global eigenvalue 
problem characterizing these latter two types of instabilities is formulated in i i4.2l 

For the equilibrium problem p.8p . it is analytically intractable to determine all possible equilibrium solution 
branches for fc-spot patterns in an arbitrary two-dimensional domain as the parameters A and D are varied. However, 
some partial analytical results are obtained in [j6] for the special case of k spots equally spaced on a circular ring 
that lies within, and is concentric with, a circular disk domain. For this special case, the Green's matrix in p.lOp is 
circulant, and the equilibrium problem is reduced to determining the equilibrium ring radius for the pattern. 

For the related Schnakenburg model, it was shown in §2.4 of |38) that near the spot self-replication threshold, 
i.e. for Sj near S2, the direction at which the spot splits is always perpendicular to the direction of the motion 
of the spot. This result was derived in [3 8) from a center- manifold type calculation involving the four dimensional 
eigenspace associated with the two independent translation modes and the two independent directions of splitting. 
Since this calculation in [38 involves only the inner region near an individual spot, it also applies directly to the GS 
model (jl.ip . This qualitative result is stated as follows: 

Principal Result 3.2: Consider the GS model (jl.ip with e <^ 1, A = 0{~elne), and t <C 0{e~'^). Suppose that 
Sj > Y.2, with Sj — S2 ^ 0^ for some unique index J in the set j ~ 1, . . . , fc. Then, the direction of splitting of the 
J^^ spot is perpendicular to the direction of its motion. 

4 Fast Instabilities of the Quasi-Equilibrium Spot Pattern 

In this section we study the stability of a /c-spot quasi-equilibrium pattern to either competition, oscillatory, or 
self-replication, instabilities that can occur on a fast 0(1) time-scale relative to the slow motion, of speed O(e^), of 
the spot locations. The stability analysis below with regards to spot self-replication is similar to that done in |38j 
for the Schnakenburg model. The formulation of a globally coupled eigenvalue problem governing competition and 
oscillatory instabilities is a new result. 

Let Ug and Ve denote the quasi-equilibrium solution of Principal Result 2.1. We introduce the perturbation 

u(x, t) = Ue + e^*ri{x) , ^(x, t) = + e'^ *0(x) , 

for a fixed spatial configuration Xi, . . . , x^ of spots. Then, from (II. ip . we obtain the eigenvalue problem 

e^A<f)- {1 + \)(t) + 2AueVe(t) + Av^rj = , L'Ary - (1 + tA)?? - 2ueUe0 - -We?] = . (4.1) 

In the inner region near the j^^ spot, we recall from (|2.ip that ^ a^/d ^^ ^^'^ ~ '^here Uj,Vj is the 
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radially symmetric solution of the core problem (|2.3p . Next, we define 



y = e-i(x-x,), v = j^N,, <l'=^^o^ (4-2) 



so that (|4.ip transforms to 

Ay$^- - (1 + A)$, + 2U,V,^, + VfN, - , AyTV, - V^N, - 2U,V,<P, = ^ (1 + rX) N, . (4.3) 

Then, assuming that D — 0(1) and r <C ©(e^^), we can neglect the right-hand side of the equation for Nj in (14.31) . 

Next, we look for angular perturbations of the form $j = e""^<&j(p) , Nj = e^™^Nj{p), where m > is a non- 
negative integer, 6 — arg(y), and p = |y|. Then, from (|4.3p . Nj{p) and $j(p) satisfy 

$"-h-$j--j^$j-(l + A)$j-f 2;7jV,$j+y/iVj=0, 0<p<c5o, (4.4 a) 
1 - . 

+ -TVj - —AT^- - y^2^^. _ 2UjVj^j = , < p < oo , (4.4 b) 

with boundary conditions 

l>^(0)==0, A^(0)=0, $j(p)->0, as p^oo. (4.4 c) 

Since the far-field behavior of Nj is different for m = and m > 2, we will consider these two different cases separately 
below. For m = 1, which corresponds to translation invariance, it follows trivially that A = is an eigenvalue of the 
local eigenvalue problem. For this translation eigenvalue, a higher-order analysis would show that A = 0{e'^) when 
£ — > 0. Since any weak instability of this type should be refiected by the properties of the Hessian of the DAE system 
of Principal Result 3.1 for the slow spot dynamics, the mode m = 1 is not considered here. 



4.1 Non- Radially Symmetric Local Perturbations: Spot Self- Replication Instabilities 

An instability of (|4.4p for the mode m = 2 is associated with the initiation of a peanut-splitting instability. Instabilities 

for the higher modes m > 3 suggest the possibility of the initiation of more spatially intricate spot self-replication 

events. Thus, the eigenvalue problem (14. 4p with angular modes to > 2 initiate angular deformations of the spot 

profile. For this range of to, the linear operator for Nj in (j4.4 allows for algebraic decay of Nj as p -H- oo owing to 

the rn^Njl p^ term. As such, for to > 2 we impose the far-field boundary condition that TV^ — as p — ?► oo . 

For TO > 2, the eigenvalue problem ()4.4p is coupled to the core problem p.3p for Uj and Vj, and can only be 

solved numerically. To do so, we first solve the EVP (|2.3p numerically by using COLSYS (cf. j2j). Then, we discretize 

(|4.4p by a centered difference scheme to obtain a matrix eigenvalue problem. By using the linear algebra package 

LAPACK [l] to compute the spectrum of this matrix eigenvalue problem, we estimate the eigenvalue Aq of (14. 4p 

with the largest real part as a function of the source strength Sj for different angular modes to > 2. The instability 

threshold occurs when Rc(Ao) = 0. We find numerically that Ao is real when Sj is large enough. In the left subfigure 

of Fig. El we plot Re(Ao) as a function of the source strength 5*^ for to = 2, 3, 4. Our computational results show that 

the instability threshold for the modes to > 2 occurs at Sj — S^, where S2 ~ 4.31, E3 ~ 5.44, and E4 w 6.14. In 

the right subfigure of Fig. [51 we plot the eigenfunction {^j,N-j) corresponding to Aq = with to = 2 at S'j = E2. 

We emphasize that since TVj — as p — ?> cxd, the initiation of a spot self- replication instability is determined 

th 

through a local stability analysis near the j spot. We summarize the result in the following statement: 
Principal Result 4.1: Consider the GS model (|l.ip with e <C 1. A = ©(— elne), and r ^ 0{e~'^). We define 
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(a) Re(Ao) vs. Sj (b) ($,,iVj) vs. p 

Figure 6. Numerical results for the principal eigenvalue Aq of (j4.4l) with mode m > 2. (a) Re{\o) vs. Sj; heavy solid 
curve is for m ~ 2 with E2 = 4.31, the solid curve is for m ~ 3 with E3 = 5.44, and the dashed curve is for m — A 
with E4 = 6.14. (b) For m = 2, the eigenf unctions {^j{p), Nj{p)) near Aq = with Sj = E2 ~ 4.31 are shown. The 
solid curve is ^j{p), and the dashed curve is Nj{p). In this subfigure the maximum value of has been scaled to 
unity. 



V and A as in (j2.5l) . In terms of A, D, and v, we calculate Si,...,Sk for a k-spot quasi- equilibrium pattern from 
the nonlinear algebraic system ()2.7|) . Then, if Sj < S2 ~ 4.31, the j^^ spot is linearly stable to a spot deformation 
instability for modes m > 2. Alternatively, for Sj > E2, it is linearly unstable to the peanut- splitting mode m — 2. 

We now show numerically that the peanut-splitting linear instability leads to a nonlinear spot self-replication event. 
This suggests that the bifurcation as Sj increases above S2 is subcritical. To show this, we formulate a time-dependent 
inner, or core, problem near a single spot, defined in terms of the local inner variables 



. Uiy,t), v^l£viy,t), y = e-i(x-x,). (4.5) 

Ay JJ c 

Then, from (jl.ip . we obtain to leading order that U and V satisfy the time-dependent parabolic-elliptic problem 
Vt ^ AyV ~V + UV'^ , AyU-UV^^O, yeR^ V^O, U^Slii\y\, as |y| ^ cx) . (4.6) 



From our eigenvalue computations, based on ()4.4p . the radially symmetric equilibrium solution to ()4.6p for U and 
V exhibits a peanut-splitting linear instability when S* > E2 ~ 4.31. To determine whether this linear instability 
leads to a nonlinear spot self-replication event when S" > E2 , we use FlexPDE (cf . 26 ) to compute solutions to (|4.6p 
in a large disk of radius |y| = Rm = 30, and with initial data 

V(y.0).|,„c,r,|y|/2,, a(y,0).l-i2:M|^^ ,„) 

The asymptotic boundary condition d\y\U = S/\y\ is imposed at |y| = Rm — 30. We set S — 4.5 > E2, and in Fig. [7] 
we plot the solution V ioT t < 300 showing a nonlinear spot self-replication event. Owing to the rotational symmetry 
of this problem, the direction of spot-splitting observed in Fig.[7]is likely due to small numerical errors or grid effects. 
In contrast, if we choose S* = 4.1 and the same initial condition, then there is no spot self-replication for ()4.6p (not 
shown). Therefore, this numerical evidence supports the conjecture that the peanut-splitting instability associated 
with the m = 2 mode initiates a nonlinear spot self- replication event when 5 > S2. 

Next, we numerically investigate spot-splitting for values of Sj that well-exceed the threshold S2 « 4.31 for 
the peanut-splitting instability. From Fig. 6(a)[ the threshold values of Sj for higher splittings are E3 ~ 5.44 and 
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(a) t = 



(b) t = 100 



(c) t = 130 



□ 



(d) t = 140 




(e) t = 170 



□ 



(f ) t = 300 



Figure 7. In a circular domain with radius Rm = 30, we compute numerical solutions to ()4.6|) by FlexPDE (cf. \2Ql ) 
using the initial condition in (|4.7p . The solution V with S — 4.5 is plotted at t — 0, 10, 100, 130, 140, 170, 300. 



6.14, corresponding to m = 3 and m = 4, respectively. From Fig. 6(a) we observe that the growth rates 



associated with these further unstable modes are comparable to that for the mode m = 2 when Sj « 6.5. 
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Figure 8. One-spot pattern in the unit square [0,1] x [0,1]. Let D = 1.0, e = 0.02 and xi = (0.5,0.7). We set 
A = 9.9213 (top row), A = 12.643 (middle row), and A = 13.285 (bottom row), so that Si « 4.5 > T,2, Si « 6.0 > S3 
and Si ~ 6.4 > E4, respectively. The numerical solutions for v from (jl.ip are computed using VLUGR (cf. 'fil), and 
are plotted at different instants in time in the zoomed spatial region [0.325,0.675] x [0.45,0.8]. 

For illustration, we consider a one-spot pattern with D = 1 and e = 0.02 in the unit square [0, 1] x [0, 1] with 
a spot centered at xi = (0.5,0.7). From (|2.21l) . we then calculate A from A = Si + 2nvRi i Si + vx{Si), where 
Ri^i is the regular part of the reduced-wave Green's function that can be calculated from (jA.ll|) of Appendix A. We 
numerically study the spot-splitting process by using VLUGR (cf. [6]) to compute solutions to (|l.ll) for A = 9.9213, 
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A = 12.643, and A = 13.285, corresponding to Si 4.5 > i;2. Si = 6.0 > S3, and Si = 6.4 > S4, respectively. 
The results are shown in Fig. |8] in the sub-region [0.325,0.675] x [0.45,0.8]. These computations show that spot 
self- replication, leading to the creation of two distinct spots, is a robust phenomena for (jl.ip whenever > S2- 
Although for Si = 6.0 and 5*1 = 6.4 the initial instability leads to a crescent pattern for the volcano profile for V , 
eventually two spots are created from this instability. Therefore, these numerical results support the conjecture that 
the unstable mode to = 2 dominates any of the other unstable modes with to > 2 in the weakly nonlinear regime, 
and eventually leads to the creation of two spots from a single spot when > S2 • 



4.2 Radially Symmetric Local Perturbations: Competition and Oscillatory Instabilities 

In §4.11 the stability of the spot profile to locally non-radially symmetric perturbations was studied numerically. In 
this subsection, we examine the stability of the spot profile to locally radially symmetric perturbations of the form 
Nj = Nj{p) and $j — ^j{p), with p = \y\, which characterize instabilities in the amplitudes of the spots. With the 
assumption that rA <C 0(e^^Z3), (|4.3p reduces to the following radially symmetric eigenvalue problem on < p < 00: 



+ -/i - + ^U,V,^, + V^'N, = X<P, , iV; + ^N^ - V^N, - 2U,V,-^, = , (4.8 a) 

(0) = N'j (0) = ; $j (p) ^ , Nj (p) ^ Cj In p + Bj + o{l), as p^oo. (4.8 h) 

Here Uj and Vj is the radially symmetric solution of the core problem (|2.3|) for the j^^ spot. From the divergence 
theorem, the constant Cj in ()4.8 b\ is given by Cj = {2UjVj^j + V^Nj) pdp. We emphasize that the operator in 
(|48a|) for Nj reduces to A^j' + p^^iVj w for p > 1, and so we cannot impose that iVj — > as p — > 00. Instead, we 
must allow for the possibility of a logarithmic growth at infinity for Nj , as written in (|4.8 &P . This growth condition, 
which will lead to a global coupling of the k local eigenvalue problems, is in contrast to the decay condition as p — > 00 
used in §4. II for the stability analysis with respect to locally non-radially symmetric perturbations. 

To formulate our eigenvalue problem we must match the far-field logarithmic growth of Nj with a global outer 
solution for rj. This matching globally couples the local eigenvalue problems near each spot. To determine the 
problem for the outer solution for 77, we use the fact that v is localized near Xj for j ~ 1, . . . ,k. Then, from (|2.ip 
and (|4.2p . we represent the last two terms in the rj equation of (|4.1|) in the sense of distributions to obtain that 
2uv(j) + v^?] - 2neVDA~^ Y^^^ Cj 5(x-Xj), where C, = {2UjVj'^>j + Vj^Nj) pdp and (5(x-Xj) is the Dirac delta 
function. Therefore, in the outer region, we obtain from (14. ip that rj satisfies 

k 

At?-^ -^^-^77 = ^Q(5(x-xj), xen- dnv = o, xedn. (4.9) 

This outer solution can be represented in terms of a A-dependent Green's function as 

^ = -:|^E^^-Ga(x;x,), (4.10) 

where G\ (x; Xj ) satisfies 

AGa- -^^-^^Ga = -(5(x-x^), xen; dnCx^O, x £ dn , (4.11a) 

Ga(x; Xj) ~ — — !- In [x — Xj[ -f i?Ajj -I- o(l) as x — ^ x^- . (4.116) 
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We remark that the regular part Rxj,j of G\ depends on Xj, D, and tA. 

The matching condition between the outer solution (|4.10l) for 77 as x — ?> Xj and the far-field behavior (j4.8 &P as 
/5 — )• 00 of the inner solution Nj near the spot, defined in terms of 77 by (|4.2p . yields that 



where Gaij = 



GA(x,j;xj) and ly 



k 



2-KE 



Cj In |x 



-1/lne. This matching condition provides the k equations 



(4.12) 



Cj (1 + + vBj + 2TTv^CiGxij = , 



j l,...,fc. 



(4.13) 



We remark that the constants r and D appear in the operator of the A-dependent Green's function defined by (|4.1ip . 

For D — 0(1) and r — 0(1), we consider the leading-order theory in v based on assuming that v = — l/lne<C 1. 
Then, to leading-order in (|4.13p yields that Cj = for j = 1, • • • , fc. This implies that Nj in (|4.8p is bounded as 
p 00, and so we can impose that Nj{p) as p 00. Therefore, when r = 0{l) and D = 0(1), then, to leading 
order in v, the eigenvalue problems (j4.8|) for j — 1, . . . ,fc are coupled together only through the determination of 
the source strengths Si, . . . , Sk from the nonlinear algebraic system (|2.7p . For this leading-order theory, we compute 
the real part of the principal eigenvalue Aq of (|4.8p as a function of 5*^, subject to the condition that Nj{p) as 
p — )■ 00. This computation is done by discretizing (|4.8p by finite difii'erences and then calculating the spectrum of 
the resulting matrix eigenvalue problem using LAPACK (cf. |lj). The plot of Re(Ao) versus Sj is shown in Fig. 9(a) 



for the range Sj < 7.5, which includes the value Sj = S2 « 4.31 corresponding to the spot self-replication threshold 
of i j4.1l In Fig. |9(b)] and Fig. 9(c)[ respectively, we plot the eigenfunctions $j(p) and Nj{p) for two different values 
of Sj. For this leading-order-in-;^ local eigenvalue problem, our computations show that Re(Ao) < for Sj < 7.5. 
Therefore, we conclude that an instability can only be generated through the global coupling of the local eigenvalue 
problems. This coupling occurs when we do not make the v ^ 1 approximation in ()4.13p . 
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(a) Re(Ao) vs. Sj 



(b) vs. p 



(c) TVj vs. p 
vs. Sj subject to the condition that Nj 



Figure 9. Left figure: The real part of the largest eigenvalue Re{Xo) of 
is bounded as p ^ 00. Middle figure: $j vs. p for Sj = 3.0 (solid curve) and Sj — 6.0 (dashed curve). Right figure: 
Nj vs. p for Sj = 3.0 (solid curve) and Sj — 6.0 (dashed curve). 



Since v = — 1/ In e decreases only very slowly as e decreases, the leading-order approximation Cj = for j = 1, . . . , k 
to (|4.13p . which does not lead to any instabilities, is not expected to be very accurate. Consequently, we must examine 
the effect of the coupling in (|4.13l) . Since (14. 8p is a linear homogeneous problem, we introduce Bj by Bj — BjCj, and 
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we define the new variables $j and Nj by $j — Cj^j and Nj = CjNj. Tlien, (I4.8P on < p < oo becomes 



1 



$;-(0) = 7Vj(0) = 0; <i'j(p)^0, 7V,(p) 



p ' 

lnp + 



2C/,y,<i>, 



0, 



as p — >■ oo . 



(4.14 a) 
(4.14 6) 



The constant -Bj in (j4.14 &P is a function of the, as yet unknown, complex-valued eigenvalue parameter A. It also 
depends on the source strength Sj through the solution Uj and Vj to the core problem (|2.3p . 
In terms of 13 j , (I4.13P transforms to the homogeneous linear system for Ci , . . . , Cfc given by 

k 

CJ{l + 2^:vR^.J^J+vB.J) + 2^vY,C^GxJ,^^Q, j = l,...,fc. (4.15) 

It is convenient to express (j4.15p in matrix form as 

Mc = 0, M=I + vB + 2TTvgx, (4.16 a) 

.,Ck)'^. In (|4.16 a\\ . ;B is a diagonal matrix and Gx is the 



where / is the k x k identity matrix and c — (Ci,.. 
A-dependent symmetric Green's matrix defined by 



I Rxi.i 

G\2A 



G\ 1,2 
R\2.2 



G\ i.k 

G\2,k 



\ G 



XkA 



G 



A ka 



G 



Xk.k 



J 



B = 



( B, 


V 





B2, 



\ 





B 



(4.16 6) 



k J 



Since G\ij ^ G\j^i when A is complex-valued, where the overbar denotes complex conjugate, Gx is not Hermitian. 

The eigenvalue A is determined from the condition that det(A^) = 0, so that there is a nontrivial solution c 7^ 
to (|4.16 ap . This condition leads, effectively, to a transcendental equation for A. The roots of this equation determine 
the discrete eigenvalues governing the linear stability of the /c-spot quasi-equilibria to locally radially symmetric 
perturbations near each spot. We summarize our formulation of the global eigenvalue problem as follows: 

Principal Result 4.2: Consider a k—spot quasi- equilibrium solution to the CS model p.ip . For e — > 0, with 
D = 0(1), A = ©(—sine), and tX ^ 0{e^'^), the stability of this pattern to locally radially symmetric perturbations 
near each spot is determined by the condition det{M) = 0, where Ai is defined in (j4.16p . The diagonal matrix B in 
(|4.16p is determined in terms of Sj and A by the local problems (I4.14p for j = 1, . . . , fc. // the principal eigenvalue 
Aq of this global eigenvalue problem is such that i?e(Ao) < 0, then the k—spot quasi- equilibrium solution is linearly 
stable to locally radially symmetric perturbations near each spot, and it is linearly unstable if Re{Xo) > 0. 

We remark that (|4.16p couples the local spot solutions in two distinct ways. First, the A-dependent terms R\jj 
and G\ij in the Green's matrix G\ in (|4.16p depend on D, on rA, and on the spatial configuration xi, . . . ,Xfc of 
spots, as well as the geometry of fl. Secondly, the constant Bj in the matrix B depends on A and on Sj. Recall that 
the source strengths Si, . . . ,Sk are coupled through the nonlinear algebraic system (j2.7p . which involves A, D, the 
reduced-wave Green's function, and the spatial configuration of spot locations. 

From our numerical study of this global eigenvalue problem in ^JSHU there are two mechanisms through which a k- 
spot quasi-equilibrium pattern can lose stability. Firstly, for k> 1 there can be a complex conjugate pair of eigenvalues 
that crosses into the unstable half-plane Re(Ao) > 0. This instability as a result of a Hopf bifurcation initiates an 
oscillatory profile instability, whereby the spot amplitudes undergo temporal oscillations. Such an instability typically 
occurs if r is sufficiently large. Alternatively, for fc > 2, the principal eigenvalue Aq can be real and enter the unstable 
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right half-plane Re(Ao) > along the real axis Im(Ao) = 0. This instability, due to the creation of a positive 
real eigenvalue, gives rise to an unstable sign-fluctuating perturbation of the spot amplitudes and it initiates a 
competition instability, leading to spot annihilation events. This instability can be triggered if D is sufficiently large 
or, equivalently, if the inter-spot separation is too small. 

The global eigenvalue problem leading to the stability formulation in Principal Result 4.2 is a new result, and 
essentially can be viewed as an extended NLEP theory that accounts for all terms in powers of v. In Appendix B, we 
summarize the main results for the leading-order-in-i^ NLEP stability theory of !62 based on the parameter range 
D = 0{i^-^) and A ^ 0{e[- \ne]^^^), and we show how our global eigenvalue problem can be asymptotically reduced 
to the leading order NLEP problem of [62] in this parameter regime for D and A. 

4.3 Symmetric Spot Patterns and a Circulant Matrix 

The global eigenvalue problem is rather challenging to investigate in full generality owing to the complexity of the 
two different coupling mechanisms in (j4.16l) . However, for the special case where the spot configuration Xi, . . . ,Xfc 
is such that Q, and consequently Qx, are circulant matrices, then the complexity of this eigenvalue problem reduces 
considerably. For this special arrangement of spot locations, the spots have a common source strength Sc — Sj for 
j = 1, . . . , k, where Sc satisfies the nonlinear algebraic equation (j2.22p . Hence, the inner problem (|4.14p is the same 
for each spot, which enforces that I3j = I3c for j = 1, . . . , k, where I3c = Bc{X, Sc)- Therefore, we can write B = I3cl 
in (|4.16p . Moreover, let v be an eigenvector of the A-dependent Green's matrix Qx with eigenvalue uj\ — ll!\{tX), 
i.e. Qxv ~ WAV. We note that the matrix Q\ is also circulant when Q is circulant. Then, the condition that M. in 
(I4.16P is a singular matrix reduces to k transcendental equations in A. We summarize the result as follows: 

Principal Result 4.3: Under the conditions of Principal Result 4--2, suppose that the spot configuration Xi, . . . ,Xfe 
is such that Q, and consequently Q\, are circulant matrices. Then, the eigenvalue condition det{A4) = for (|4.16p 
reduces to the k transcendental equations for A given by 

fj = l + iyBc + 2ni'ujxj{TX) =0, (4.17) 

where uj\j (tA) for j — 1, . . . , k is any eigenvalue of the matrix Q\ . The k-distinct eigenvectors v of Qx determine the 
choices for c = (Ci, . . . , Ck)^ . By equating real and imaginary parts, J^.JTp can he reduced to 

Re{fj) = l/i/ + Re(Bc) + 27rRe(wAj(TA)) = , Im(/j) = Im(i?c) + 27rIm(wAj(rA)) = . (4.18) 

When Qx is a circulant and symmetric matrix, then its spectrum, as needed in Principal Result 4.3, can be 
determined analytically. To do so, let the row vector a = (oi, . . . , a^) denote the first row of Q. Since Qx is circulant it 
follows that all of the other rows of G can be obtained by cycling the components of the vector a. In addition, since Gx 
is also necessarily a symmetric matrix it follows that 02 = a^, 03 — afc_i, . . ., and Oj — ak+2-j for j — 2, . . . , \k/2~\ , 
where the ceiling function \x~\ is the smallest integer not less than x. We recall that if a fc x fc matrix is circulant, its 

th 

eigenvectors Vj and eigenvalues ujxj, which consist of the k roots of unity, are given by 

cux, = , V, = (1, e^-O-DA , . . . , e^M,-i)(k-i)/k-^T j ^ 1, . . . , fc . (4.19) 

m=0 

Here Um is the component of the row vector a. 
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For illustration, let fc = 3. Then, the eigenpairs are vi = (1, 1, 1)-^ with lu\i = ai +02 + 03, V2 = (1, e^'^*/'^, e^'^'/-^)^ 
with WA2 = ai +0262^^*/=^ + a3e''^*/3^ and V3 = (1, e^""''^, e^-^'^^f with a;A3 = ai + a2e^'^'^^ + ase^''''/^ . Then, since the 
matrix is also symmetric, we have 02 = 03, so that uj\2 — = 01 — 02, which yields one eigenvalue of multiplicity 
two. Then, since any linear combination of V2 and V3 is also an eigenvector, we take the real part of V2 as one such 
vector, and the imaginary part of V2 as the other. In summary, we can take Vi — (1, 1, 1)-'", V2 — (1, —0.5, —0.5)"^ 
and V3 = (0, V3/2, as the three eigenvectors of Qx for the three-spot pattern. 

In general, the symmetry of Qx implies that aj — ak+2-j and Vj — 'Vk+2-j, and Vjm — ''^j(fe+2-m) for 'tt, = 
2, . . . , [fc/2] and j — 2, . . . , [fc/2] . Here Vjm denotes the rrfi'- component of the column vector vj. Since the eigenvalue 
is LUxj — Vj - a, we have ujxj ~ ^\{k+2-j) ■ This generates [fc/2] — 1 eigenvalues with multiplicity two, whose eigenvectors 
can be obtained by taking any linear combination of two complex conjugate eigenvectors. For instance, for the 
eigenvalue loxj the corresponding two eigenvectors can be taken as the real and imaginary parts of Vj , respectively. 
In summary, we know that all of the eigenvectors can be chosen to be real, but that the eigenvalues in general will be 
complex when A is complex. This leads to the following result for the spectrum of the fc x fc symmetric and circulant 
Green's matrix Qx whose first row vector is a = (oi, . . . , a™): 



uJxj 



— J2m=o '^o^ f^lU—lllIl \ a„i+i, eigenvalues with multiplicity 2, 
= (l,c"os(H^),...,cos(^.(^K^)), ^'-''^ 



vr,_. =fo,sin ,...,sin ^^(l^lfcii) ) , 2, . . . , rfc/2l + 1 . 



k+2-j k ;>--->^^^^\^ k 

Note that if fc is even, then a;A( [fc/2] +1) = X)m=i(^f )'"^^^™ a simple eigenvalue with eigenvector (1, —1, • • • ,1, — 1)'^. 

For a symmetric spot pattern, we can simply substitute ()4.20p into (|4.18p of Principal Result 4.3 to derive the 
transcendental equation associated with the eigenvector Vj of Qx- This yields, 

fc k 

u-^ +Re{B^) + 2Tr ^Re(am)vjm^O, Im(4) + 27r ^ lm(a„0^'jm = , j = 1, . . . , rfc/2] + 1 . (4.21) 

m—l m—1 

Here the row vector a has components ai — Rx 1,1, and a„j = Gx i,m for m = 2, . . . , fc, where Gx i,m and Rx 1,1 are 
determined from the A-dependent Green's function and its regular part satisfying (|4.1ip . In (I4.2ip . vjm denotes the 
TO^^ component of the column vector Vj. We remark that when A is real, (|4.2ip can be reduced to only one equation. 

At the onset of a competition instability, where an unstable eigenvalue first enters the unstable right-half plane 
Re(A) > along the real axis Im(A) = 0, the result (I4.2ip can be simplified further. Upon differentiating the problem 
for Uj and Vj in (j2.3p with respect to Sj, and then comparing the resulting system with the eigenvalue problem in 
(|4.14p , it readily follows that 

h = ^V,, N.^^U,, when A = 0. (4.22) 

Therefore, for the special case where 5j = 5c for j = 1, . . . , fc, the constant Be in (|4.21l) when A = can be calculated 
in terms of the derivative of x by 

B^^x'iSc), when A = 0. (4.23) 



A plot of x'i^c) is shown below in Fig. 10(a) Thus, for a configuration of spots for which the Green's matrix is 
circulant symmetric, the threshold condition for a competition instability, corresponding to setting A = in (|4.2ip 
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and recalling (|2.22l) . is to solve the coupled system 

x\Sc) + 2Trujxj ^ -u-^ A = Sc{l + 2Tny9) + i^x{Sc); A^e-^i^AVd, = ]^ ■ (4-24) 

Since G and Q\ coincide at A = 0, then in (|4.24p is the eigenvalue of G with eigenvector e = (1, . . . , 1)-^, i.e. Qe — 9e, 
and u}\j for j = 1, . . . , [fc/2] + 1 is any of the other eigenvalues of G as given in (|4.20p when A = 0. We remark that 
when k is even, our computational results in ^JS] - [H] below will show that the most unstable mode for a competition 
instability is the sign-fluctuating mode (1, —1, • • • ,1, —1)^, which corresponds to setting j — \k/2] + 1 in (I4.24p and 
using ujx{[k/2]+i) = STO=i(^l)™"^"m- Here ai = and a™ = Gi,™ for to = 2, ... , k, is the first row of G- 

Owing to the considerable reduction in complexity of the global eigenvalue problem when G and G\ are circu- 
lant matrices, in ^JSHS] we will only compute competition and oscillatory stabilities thresholds of quasi-equilibrium 
spot patterns that lead to this special matrix structure. The numerical approach used to compute these thresholds 
associated with the stability formulation of Principal Result 4.3 is outlined in Appendix C. 




Finally, we remark that the numerical implementation of our asymptotic theory in fJ2] and 21 for spot dynamics 
and stability relies on the computation of the three functions xi^), x'i^), and 7(5*), defined in terms of the core 
problem by (|2.3p and (13. 7p . The result for x(5') was given previously in Fig. 1(b) [ while x'(5') and 7(5) are plotted 
in Fig. 10(a) and Fig. |10(b")| respectively. 



5 One- and Two-Spot Patterns on the Infinite Plane and in Large Domains 

In this section we apply the asymptotic theory of [J2] - 14] to the special case of either a one- or a two-spot pattern 
on the infinite plane . For this problem we can set _D = 1 in (|l.ip without loss of generality. When i7 = , the 
reduced-wave Green's function and its regular part satisfying (jl.2p with _D = 1 is simply 

G(x;x,)-^ifo(|x-x,|) , i?^.^. ^-L(ln2-7e). (5.1) 

Here 7e is Euler's constant, and KQ{r) is the modified Bessel function of the second kind of order zero. 

A fc-spot quasi-equilibrium solution with spots at Xj € for = 1, . . . , fc is given in (|2.8p of Principal Result 2.1. 
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From (|2.7p and (|5.ip , the source strengths Si, . . . , Sk satisfy the nonhnear algebraic system 

k ^ 
^ = S'j(l + i^(ln2-7,)) +z/^5,i^o|x. --Xj-|+i/x(5'j), i = l,...,fc; y^^e^VA, 1, = -—. (5.2) 

To determine the stability of the fc-spot quasi-equilibrium solution to locally radially symmetric perturbations, 
we must find the eigenvalues of the global eigenvalue problem of H4.2I consisting of (|4.14p coupled to (|4.15p . The 
A-dependent Green's function and its regular part, satisfying (|4.1ip with Z) = 1 in M'^, are 

Ga(x;Xj) = ^A^o (VTT7A|x-x,|) , = ^(ln2 - 7^ - logN/TT^) . (5.3) 

Since A can be complex, we must take log z to be the principal branch of the logarithm function, and choose for 
z = vT+tA the principal branch of the square root, in order that Kq{z) decays when z — > cx) along the positive real 
axis. In terms of this Green's function, the homogeneous linear system (|4.15p becomes 

k 

C.J (l + v{\n 2 - 7e - log Vl + rA + 4)) + ^'^^ (Vl + tA |x,; - x^ |) = , j = l,...,k. (5.4) 

Below, we consider the specific cases of one- and two-spot solutions in . The numerical procedure used below to 
compute the competition and oscillatory instability thresholds is summarized in Appendix C. 



5.1 A One-Spot Solution in 



For a one-spot solution with spot at the origin, ()5.2p reduces to the scalar nonlinear algebraic equation 

A^Si+vSi{\n2--f^) + vx{Si), A = e~^vA, v ^ . (5.5) 

me 

In order that Ci 7^ in (|5.4p . we require that A satisfy the nonlinear transcendental equation 

+ (ln2-7e-logVl + rA) = 0. (5.6) 

Here B\ ~ B\{\,S\) is to be computed from (|4.14p in terms of the solution [/i and V\ to the core problem (|2.3p . 
Our computational results lead to a phase diagram in the A versus t parameter-plane characterizing the stability of 
a one-spot solution. In the numerical results below we have fixed e — 0.02. 

From (j5.5p we compute that there is no quasi-equilibrium one-spot solution when A > Af = 0.1766. This existence 
threshold for A, which was obtained by varying Si in (|5.5p on Si G [0.22, 7.41], is shown by the lower thin solid line 



in Fig. 11(a) From ^4.1[ the spot self-replication threshold Ag ~ 0.311 for A is obtained by setting 6*1 = E2 ~ 4.31 



and x(^2) ~ —1.783 in (15. 5p . It is shown by the upper horizontal dotted line in Fig. 11(a) Moreover, in Fig. 11(a) 



the solid curve is the Hopf bifurcation threshold A versus th as computed from ()5.6p . At this Hopf bifurcation value, 
a complex complex conjugate pair of eigenvalues first enters the right half-plane. The resulting parameter-plane as 



shown in Fig. 11(a) is divided into four regions with the following solution behavior. In Regime a with A < Af, the 
quasi-equilibrium solution does not exist; in Regime /3, enclosed by the thin solid horizontal line Af = 0.1766 and 
the heavy solid Hopf bifurcation curve, the quasi-equilibrium solution is unstable to an oscillatory profile instability; 
in Regime C it is stable; and in Regime 6, it is unstable to spot self-replication. In Fig. ll(b)[ we fix ^ = 0.18 and we 
show that the real part Re(A) increases with t on the range r S [30, 40]. Note that for A = 0.18, the Hopf bifurcation 
threshold th — 36.06 is where the complex conjugate pair of eigenvalues intersect the imaginary axis. 

For three values of e, in Fig. [T^ we compare the stability threshold obtained from (|5.6p with that obtained from 
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Figure 11. One-spot solution in R'^ with e = 0.02. (a) Phase- diagram in the parameter-plane A vs. t. The solid 
curve shows the Hopf bifurcation threshold th , the lower horizontal thin solid line plots our existence threshold 
Af 0.1766, the dot-dashed horizontal line is the existence threshold Af^ 0.1758 from NLEP theory, and the 
upper heavy dotted horizontal line indicates the spot- splitting threshold As = 0.311. Regime a: No quasi- equilibrium 
solution exists; Regime 13: Oscillations in the spot amplitude; Regime Stable one-spot solution; Regime 9: spot 
self-replication regime, (b) Fix A = 0.18, the spectrum in the complex A plane is shown for t e [30,40]. At t = 30, 
A = -0.026±0.296i. AtT = 36.05, A = -0.000049±0.285 i. AtT = 40, X = 0.015±0.278 i. As r increases a complex 
conjugate pair of eigenvalues enters the right half-plane. 
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Figure 12. One-spot solution in R^. In the A vs. j = — Inr/ Ine plane, we plot our existence threshold Af by the thin 
dashed horizontal line, our Hopf bifurcation threshold A by the heavy solid curve, and the stability threshold Asw from 
the NLEP theory by the dot-dashed curve (a) e = 0.02, Af = 0.17659; Af^ = 0.17575; (b) e = 0.01, Af = 0.095875, 
Af^ = 0.095342; (c) e = 0.005, Af = 0.051432, A/^ = 0.051133. Here Af^ is the existence threshold from NLEP 
theory given in (|5.7p . 



the NLEP analysis of [60], as summarized in Appendix B. From (IB.3[) and (|B.5p of Appendix B, there are two 
threshold values Af^ and Agy^ predicted from the leading-order- in-i^ NLEP theory of 1601 . For A < Af^ a one-spot 
quasi-equilibrium solution does not exist. For A < Ag^ this solution is unstable to an oscillatory instability. These 
threshold values, with Asw depending on t and the ground-state solution w of (IB.1|) . are 
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(5.7) 



For e = 0.02, we get Afw ~ 0.1757, which agrees very well with our existence threshold Af w 0.1766. In Fig. [T2l we 
compare the NLEP stability threshold As^, versus 7 = — Inr / Ine with our Hopf bifurcation threshold A versus 7 = 
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— Inr/f/lne. We recall from Appendix B that the NLEP theory predicts instability when A < Asw for < r < £~^, 
but has no conclusion regarding stability or instability ii A > Asw Our Hopf bifurcation threshold predicts that 
stability changes when we cross the heavy solid curves in Fig. [T^l 

Next, we study the possibility of real- valued unstable eigenvalues that have crossed along the real-axis Im(A) = 
into the unstable half-plane Re(A) > 0. The threshold condition for this is to set A = and 13i — x'{Si) in (|5.6p to 
obtain that satisfies x'li^i) + 1^^^ +ln2 = 7e. By differentiating (|5.5p for A with respect to Si, we obtain that this 
threshold coincides with the minimum value 5*1 ,„ of 5i for the curve A versus Si as defined by (|5.5p . For S > Sim, 
we compute numerically that Re(A) < 0, while Re(A) > for S" < Sim- The critical value Sim, for which A = 0, 



corresponds to the fold point in the bifurcation diagram of u(0) versus A (such as shown in Fig. 5(b) 



5.2 A Two-Spot Solution in 

Next, we consider the case of two spots centered at xi — {a, 0) and X2 = (—a, 0), where a > and a ^ 0{e). We 
look for a symmetric solution where Si — S2 — Sc- Then, from ()5.2|) . Sc satisfies the nonlinear algebraic equation 



^ = [1 + J^(ln 2 - 7e) + lyKa (2a)] + ly xiSc) 



A 



-^uA. 



-1/lne. 



(5.8) 



For e = 0.02, we first discuss the existence of the quasi-equilibrium two-spot solution. In Fig. 13(a) we plot A 



versus Sc for three different values of a, showing the existence of a fold point Af in the graph of Sc versus A. This 
fold point Af{a) is shown in Fig. 13(b) by the dotted curve. The solid curve in Fig. |13(b)| is the corresponding 
two-term asymptotic result for Af obtained by expanding (|5.8I) in powers of v as in (|2.23p . This result shows that a 
quasi-equilibrium two-spot pattern exists for a spot separation distance of 2a only when A > Af{a). 




< 3 



1 2 3 4 5 6 



(a) A vs. Sc 



(b) Af 



a 

(c) Sc vs. a 



Figure 13. Two spots at (±a,0) in with e = 0.02. (a) A vs. Sc for a = 0.1 (solid curve), a — 0.5 (dashed curve) 
and a — 1.0 (heavy solid curve), (b) The fold point A f vs. a; two-term asymptotic result (solid curve), and numerical 
result from (|5.8|) (dotted curve), (c) Sc vs. a for A = 3.5 (heavy solid curve), A = 4.0 (solid curve) and A = 4.5 
(dashed curve). 



Next, we discuss the possibility of self-replicating instabilities for an inter-spot separation of 2a. In Fig. 13(c)[ we 
plot Sc versus a ioi A — 3.5 (solid curve), A — 4.0 (dashed curve) and A — 4.5 (heavy solid curve). Notice that for 
a fixed A> Af, there are two solutions for Sc for a given a. Since the solution branch with the smaller value for Sc 



is unstable, we only plot the large solution branch for Sc in this figure. Fig. 13(c) shows that for a fixed a, Sc is an 



increasing function of A, whereas for a fixed A, then Sc increases as a increases. For the solid curve in Fig. 13(c) 
with A — 3.5, Sc < S2 for all a, so that the spots never split for any inter-separation distance. Alternatively, the 



dashed curve for A — 4.0 and the heavy solid curve for A — 4.5 in Fig. 13(c) intersect the spot-splitting threshold 
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5c = S2 at a « 1.81 and a « 0.46, respectively. This threshold initiates a spot-replication event. The self-replication 
threshold As versus a is obtained by setting Sc — 4.31 in (j5.8|) . Moreover, from (|5.8|) . we readily observe that As 
is a decreasing function of a and asymptotes to its minimum value Asm = 3.98 as a — > 00. This minimum value is 
obtained by setting Sc = S2 and Ko{2a) — in (|5.8p . We conclude that for any A > Asm, a two-spot pattern on 
will be linearly unstable to a peanut-splitting instability if the inter-spot separation distance 2a is sufficiently large. 
Thus, for A > Asm, a spot self- replication instability is, essentially, under-crowding instability, that is triggered only 
when the two spots are too isolated from each other. In terms of the original feed-rate parameter A, the spot-splitting 
threshold As = eAs/v is given by the heavy solid curve in the phase diagram of Fig. 14(a) 



Next, we compute the stability thresholds for either competition and oscillatory instabilities. The A-dependent 
Green's function is given in ()5.3|) . Since the A-dependent Green's matrix Q\ is circulant, then Principal Result 4.3 
applies. The eigenvectors and eigenvalues uj\j of Q\ for i = j, 2 are 

1 



WA 1 = uj\i{t\) = Gxi2 + -Rah = 

WA2 = 0J\2{tX) = Rxil — G\12 = 



27r 
1 

2^ 



(In 2 - 7e - log V1 + tA) + Ko{2a^/T+T\) 



vi^(l,ir, 
V2^(l,-ir. 



(In 2 - 7e - log VTTtA) - Ko{2a\/l + TX) 

Then, from (I4.17P of Principal Result 4.3, we must determine the roots A of the two transcendental equations 

j/-i+4 + 27rwAj(rA) = 0, j = l,2, (5.9) 

where -Bc(A, Sc) is to the common value for Bj for j ~ 1,2 in (j4.14|) . 

Our numerical results below show that in some region of the A versus a parameter-plane, the two-spot solution 
becomes unstable when r is increased due to the creation of an unstable eigenvalue A from uj\i as a result of a Hopf 
bifurcation. Since this oscillatory instability is associated with the eigenvector vi = (1,1)* — {Ci,C2)^ , it leads to 
the initiation of a simultaneous in-phase oscillation in the amplitudes of the two spots. In another region of the A 
versus a parameter plane we will show numerically that an unstable real positive eigenvalue A for a;A2 can occur, 
even when t = 0, if the inter-spot distance 2a is below some threshold. Since this instability is associated with the 
eigenvector V2 = (1, —1)* — (Ci, G2)'^ , it leads to the initiation of a sign-fluctuating instability in the amplitudes of 
the two spots, which leads ultimately to the annihilation of one of the two spots. Since this instability is triggered 
when 2a is sufficiently small, it is also referred to as a competition or an overcrowding instability. The threshold 
for a competition instability associated with the dominant sign-fluctuating mode V2 = (1,-1)^ is obtained from 
the formulation (|4.24p . The threshold condition for this instability is reduced to finding the curve A = A{a), by 
eliminating the source strength Sc from the following nonlinear algebraic system: 

e 



x'(5c)+ln2-7e-Xo(2a) = -i.- 



A 



{Sc [1 -t- (ln2 - 7e + A'o(2a))] + vx{Sc)) ■ 



(5.10) 



We fix e = 0.02 in the computations below. In Fig. 14(a) the two-spot existence threshold A/ as a function of 
a e [0.02, 2.02] is shown by the lower solid curve, while the spot self-replication threshold As versus a discussed 
above is shown by the top heavy solid curve. To obtain Af we determined the minimum value of A as Sc is varied 
in ()5.8|) . In Fig. 14(a) the middle dotted curve is the threshold for the onset of a competition instability, obtained 



from (|5.10l) . These three curves separate the A versus a parameter-plane of Fig. 14(a) into four distinct regions. 
In Regime tr, the quasi-equilibrium two-spot solution does not exist. In Regime /3, the two-spot quasi-equilibrium 
solution undergoes a competition instability for any r > 0. This instability ultimately leads to the annihilation of 
one of the spots. We note that the competition threshold approaches the existence threshold when a becomes large. 
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This is because when a 3> 1 the spots are uncoupled, and there is no competition instabihty for a one-spot solution 
in R^. In Regime C, the two-spot pattern is stable to a competition instability, but becomes unstable to an oscillatory 
profile instability when r is large enough. As t increases above a certain threshold r//, a Hopf bifurcation occurs and 
the two-spot solution is unstable to the dominant in-phase, or synchronous, spot amplitude oscillation. In Regime 9 
the two-spot quasi-equilibrium solution is unstable to spot self-replication for any t > 0. 
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Figure 14. Two-spot solution m M^. (a) We plot the existence threshold Af vs. a (lower curve), the spot- splitting 
threshold As vs. a (top curve), and the competition instability threshold (dotted curve), (b) Hopf bifurcation threshold 
thi for in-phase oscillations (bottom of each pair of curves) and th2 for out-of-phase oscillations (top of each pair of 
curves), are plotted versus a for A = 0.18 (bottom pair), A — 0.2 (middle pair), and A = 0.22 (top pair). The dotted 
portions on t^i correspond to the parameter regime /3 in suhfigure a) where competition instabilities occur for any 
T. (c) Imaginary part Xi of X vs. a at the Hopf bifurcations thresholds for A = 0.18 (bottom), A = 0.2 (middle), and 
A — 0.22 (top). The heavy solid curves are the eigenvalues for in-phase oscillations, and the thin solid curves are the 
eigenvalues for out-of-phase oscillations. The dotted portions lie on thi where a competition instability occurs, (d) 
Eigenvalue path in the complex plane for A = 0.18 and a — 1 on the range r G [10, 20], with the arrow indicating the 
direction as r increases. The eigenvalues associated with in-phase (solid curves) and out-of-phase (dashed curves) 
oscillations enter the right half plane at tri and th2 , respectively. 



For three values of A, in Fig. |14(b)| we plot the Hopf bifurcation thresholds thi and th2 associated with the 
eigenvectors vi = (1,1)^ and V2 = (1,-1)^, respectively. From this figure we observe that in the limit a — > oo, 
corresponding to two isolated spots, the thresholds thi and th2 approximate the Hopf bifurcation threshold th of 
an isolated one-spot solution. 

From Fig. |14(b)[ the Hopf bifurcation curves for thi , which each have a dotted portion, end at the values a/ = 
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0.888, 0.206, 0.068 for the curves A — 0.18, 0.2, 0.22, respectively. These critical values of a correspond to the existence 
threshold A = Af{a) for the quasi-equilibrium solution. Moreover, the Hopf bifurcation threshold th2 terminates 
at another set of critical values ac = 0.904,0.28,0.16 for the curves A — 0.18,0.2,0.22, respectively. The reason for 



this disappearance is seen in Fig. 14(c) where we plot the imaginary part of the eigenvalues at the Hopf bifurcation 
thresholds for vi — (1, 1)-^ (curves with dotted portions) and V2 = (1, —1)"'" (curves without dotted portions). From 
this figure it is clear that as a —> ac from above, the complex conjugate pair of eigenvalues associated with th2 merge 
onto the real axis at the origin, which is precisely the crossing point for the onset of a competition instability. 

When the spots are too close in the sense that af < a < then there is a positive real (unstable) eigenvalue 
for any r > 0, which initiates a competition instability. If in addition, r < thi, it is the only unstable eigenvalue. 
In Fig. |14(b)] and Fig. 14(c) the dotted segments of the curves correspond to the range of a where this competition 
instability occurs. Finally, we observe from Fig. [14(b)] that THi < TH2- Therefore, for each a > ac, the in-phase 
synchronous oscillatory instability of the spot amplitudes is always the dominant instability as r is increased. 

For A = 0.18, a = 1.0, and 10 < t < 20, in Fig. |14(d)| we plot the path of the eigenvalues in the complex plane, 
corresponding to both vi (solid curves) and V2 (dotted curves), as r is increased. For each case, Re(A) increases as 
T increases. The two Hopf bifurcation values are, respectively, thi = 13.1 and th2 = 15.6. 
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Figure 15. Two-spot solution in M^. Fix A — 0.22 and e = 0.02. The Hopf bifurcation threshold thi vs. a for a 
synchronous oscillatory instability obtained from the full stability problem (|5.9p (solid curve) is compared with that 
obtained from the approximate system 15.11]) resulting from the large-T approximation (dotted curve). 




Next, we derive an approximation for the Hopf bifurcation threshold valid for t ^ 1. For r ^ 1 and a not too small, 
we exploit the exponential decay of the modified Bessel-function for large argument to obtain that the A-dependent 
Green's matrix Q\ is essentially a diagonal matrix. Therefore, in this large-r limit we can substitute waj ~ ^Aii for 
i — 1,2 into ()5.9p . In this way, and by using (|5.3|) for i?Aii, we obtain that the limiting Hopf bifurcation threshold 
th and the critical eigenvalue A = iXi are the solutions to the approximate system 

Im(4)-J, Re(4) +ln2-7e-iln(TA,) = -i^-^ (5.11) 

For £ — 0.02 and A = 0.22, in Fig. [15] we compare the threshold obtained from (|5.1ip with that obtained from the full 
transcendental equation (|5.9p (as shown initially in Fig. |14(b)| . This comparison, on the range a > 0.16 for which 
no competition instabilities occur, shows that (|5.11|) gives a good approximation provided that a is not too small. 

Finally, we derive the ODE system that determines the slow dynamics of the two-spot pattern when Sc < S2. 
Provided that no other profile instabilities are present, and that A > Af{a), we can explicitly evaluate the terms in 
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((3J)) of Principal Result 3.1 to obtain the DAE ODE system 

^--7(^c)«(2a), l{Sc)^ ^oo^l, ■ , i^eh, (5.12) 
di (j)*V^pdp, 

with Sc determined in terms of a by (|5.8p . Since K'^lr) < for any r > and [/'^□(r)! is a decreasing function of r, 
(|5.12p shows that the two spots are repelling and that their common speed ac = \ -^\ is a decreasing function of a. 
Since Ko(2a) is an exponentially small when a ^ 1, the two-spot dynamics (j5.12p is metastable in this limit, and is 
similar to the metastable repulsive dynamics studied in |25j for a two-spot solution to the GM model in R^. 

Next, we discuss the possibility of dynamically-triggered instabilities induced by the slow evolution of the two-spot 
pattern. With regards to spot self-replication, we observe that since the two spots are dynamically repelling, and Sc is 



a monotonically increasing function of a (see Fig. 13(c) I, it follows that a dynamically-triggered spot self-replication 
instability will occur when A > As7n ~ 3.98. To illustrate this let A — 4.0. Then, the dashed curve of Sc versus a 



in Fig. 13(c) shows that Sc > S2 only if a > 1.81. Now consider an initial two-spot quasi-equilibrium solution with 
a = 1 at t = 0. For this initial value of a, then A > Af (see Fig. |13(b)| . Then, since a' > and a — > 00 as i — > 00 
from the ODE (|5.12p . it follows that a — 1.81 at some long time t = 0{e^'^), which triggers the spot self-replication 
instability. From Principal Result 3.2 of ^21 the spots will split in a direction perpendicular to the x-axis. 

Thus, for the GS model (jl.ip . we conclude that there is a wide parameter range for which two spots in will 
undergo a dynamically-triggered spot self-replication event that is induced by their slow repulsive dynamics. In 
contrast, we remark that since the Schnakenburg model of ;38; involves the Neumann Green's function rather than 
the reduced-wave Green's function, the infinite-plane problem for this related RD model is ill-posed. 

Finally, we remark that the monotone increasing behavior of the Hopf bifurcation threshold thi versus a (see 
Fig. |14(b)| ), coupled to the repulsive spot dynamics of the two-spot pattern, precludes the existence of any dynamically- 
triggered oscillatory instability for a two-spot solution in R^. Namely, consider an initial two-spot quasi-equilibrium 
with initial inter-separation distance 2ao, and with a fixed r with r < THiiao), where ao > olc- Then, since a 
increases as t increases, and thx is monotone increasing in a, it follows that r < r/fi(a) for any i > 0. In a similar 
way, we can show that there are no dynamically-triggered competition instabilities in M^. 



5.3 Two-Spot Patterns in a Large Square Domain 

Next, we consider a two-spot pattern in the square domain [0, L] x [0, L\ with spots located symmetrically about the 
midline of the square at Xi^2 — {L/2±a, L/2), where a > 0. For this arrangement, Q in (|2.10p is circulant symmetric, 
with matrix entries that can be numerically calculated from (|A.11I) of Appendix A. We fix the parameter values 
D = 1, A = 4.2, and e = 0.02. Since the limiting solution for L — >■ 00 reduces to the infinite-plane problem, a large 
square domain can be used to approximate a two-spot evolution for the infinite-plane problem of ^5.2[ Our goal here 
is to examine the asymptotic prediction of spot self-replication. 

In Fig. [16] we plot Sc versus a for a square of side-lengths L — 6,8, 10. By symmetry, the equilibrium spot locations 
are at ae = L/A. In this figure, we also plot the infinite-plane result for Sc versus a, obtained from ()5.8p . The resulting 
four curves of Sc versus a are essentially indistinguishable on the range a L where the finite-boundary effects are 
insignificant. From this figure, we observe that the spot-splitting threshold Sc — '^2 ^ 4.31, represented by the thin 
solid line, intersects each of the four curves at some critical values of a, which initiates a spot self-replication event. 
We observe from Fig. [11] that Sc < S2 when a becomes close to L/2. Therefore, for a spot near an edge of the square, 
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the Neumann boundary condition effectively introduces an image spot outside of the square in close proximity to 
the interior spot, which then eliminates the possibility of spot self-replication. 




Figure 16. Comparison of two-spot patterns in a square domain and in an infinite domain. Fixing e = 0.02, D = \ 
and A — 4.2, we plot Sc vs. a for an infinite domain (heavy solid curve), and for square domains with side-lengths 
L = 10 (dashed curve), L = 8 (dotted curve), and L = 6 (dot-dashed curve). The thin solid line is the peanut- splitting 
threshold E2 ~ 4.31, which intersects the four curves at as ~ 0.775 for the infinite domain, w 0.775 for L — 10, 
as2 « 0.780 for L = 8, and a,3 ~ 0.805 for L = 6. 

For L = 8, we now compare our asymptotic prediction for spot self-replication with corresponding full numerical 
results computed from (|l.ll) using VLUGR (cf. [6]). We take an initial two-spot pattern with a — 0.4 < at time 
t — 0. At this initial time, we calculate from (|2.22|) that Sc ~ 3.77, so that the initial two-spot pattern is stable 
to spot self-replication. Then, since the equilibrium state is at Ue = 2.0 and the two-spot dynamics is repulsive, 
it follows that a increases as t increases and will eventually exceed the spot-splitting threshold a^a ~ 0.780. This 
prediction of a dynamically-triggered instability, which very closely approximates that for the infinite-plane problem 
of tj5.21 is confirmed by the full numerical results shown in Fig. 1171 




Figure 17. A two-spot pattern in the square [0,8] x [0,8]. Let A = 4.2, D = 1.0, e = 0.02 with initial spot locations 
Xi,2 (4±0.4, 4) at time t — 0. Then, Sc ~ 3.77 att — 0, but Sc > ^2 ~ 4.31 at a later time. The numerical solution 
to (jl.ip . represented as a gray-scale image of v, is plotted in the zoomed region [2, 6] x [2, 6] at times t ~ 1, 41, 131, 181, 
231,281,331,431. The spots undergo splitting in a direction perpendicular to their motion consistent with Principal 
Result 3.2. 

6 Symmetric Spot Patterns in the Square and Disk 

In this section, we implement the asymptotic theory of S|2H4]to study competition, oscillatory, and spot self-replication 
instabilities for various spot patterns for which the Green's matrix Q in (j2.9|) has a circulant matrix structure. The 
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thresholds for these instabilities are favorably compared with results from full numerical simulations obtained by 
using the PDE solver VLUGR (cf. [6]) for a square domain and the finite element code of 36 for the disk. 

In all of the numerical simulations below and in fJZIwe have set e = 0.02 so that v ~ —1/ Ine ~ 0.2556. For given GS 
parameters A and and for an initial configuration of spot locations xi, . . . ,Xfc, we compute the source strengths 
5i , . . . , S'fe from the nonlinear algebraic system (|2.7p . The initial condition for the full numerical simulations is the 
quasi-equilibrium solution (12. 8p with the values for V,(0) as plotted in Fig. l(b)[ corresponding to the computed 
values of Sj. Since this initial condition provides a decent, but not sufficiently precise, initial /c— spot pattern, we 
only begin to track the spot locations from the full numerical simulations after the completion of a short transient 
period. To numerically identify the locations of the spots at any time, we determine all local maxima of the computed 
solution V by identifying the maximal grid values. This simple procedure is done since it is expensive to interpolate 
the grid values to obtain more accurate spot locations at every time step. In this way, both the spatial trajectories 
and the amplitudes of the spots are obtained from the full numerical results. In many of the numerical experiments 
below, we show a gray-scale contour plot of the numerical solution for v. The bright (white) regions correspond to 
the spot regions where v has a large amplitude, while the dark region is where v is exponentially small. 

6.1 A Four-Spot Pattern in the Unit Square 

We first consider a four-spot pattern in the unit square VI — [0, 1] x [0, 1] with e — 0.02 fixed. The four spots are 
centered at the corners of a square symmetrically placed inside 51 as shown in Fig. 1201 For this configuration, the 
Green's matrix Q is circulant symmetric, so that there is a solution to ()2.10|) for which the spots have a common 
source strength Sc- From (|2.22[) . this common source strength Sc satisfies the nonlinear algebraic equation 

A= Sc + 2t:v6Sc + vx{Sc), Ae/{v^/D), i^ = -l/lne, 61 = + Gi,2 + Gi,3 + Gi,4 • (6.1) 

Numerical values for Rn and Gij for j = 2, . . . , 4 are obtained from the reduced-wave Green's function for the unit 
square as given in (jA.llI) of Appendix A. 

We let r denote the distance from each spot to the center (0.5, 0.5) of the square, so that r e (0, 1/V2). For D = l, 
in Fig.[TH]we plot the phase diagram of A versus r for 0.1 < r < 0.7. In this phase diagram, the spot self-replication 
threshold is obtained by setting x(^2) ~ —1.79 and S2 ~ 4.31 in ()6.1|) . The existence threshold for the four-spot 
quasi-equilibrium pattern is obtained by determining, for each fixed r, the minimum point of the curve A versus Sc 
from (|6.ip . Finally, as similar to the case of a two-spot pattern on the infinite plane as studied in H5.2[ the competition 
instability threshold in Fig. [18] is set by the sign-fluctuating eigenvector. Therefore, from (|4.24p this threshold of A 
versus r is obtained by numerically solving 

4 

x'(5c) + 2^^(-l)™-ia,„ = -z/-i; ai = i?ii ; = Gi, , j = 2,...,4, (6.2) 



together with (|6.ip . Recall that a plot of x'('S'c), as needed in (|6.2|) . was shown in Fig. 10(a) 

The phase diagram of Fig. [18] consists of four distinct parameter regimes. In Regime cr the four-spot quasi- 
equilibrium solution does not exist. In Regime /3 the quasi-equilibrium solution exists but is unstable to a competition 
instability. In Regime C the solution is unstable to an oscillatory profile instability when r exceeds a Hopf bifurcation 
threshold th- In Regime 6 the solution is unstable to spot self-replication. 

Next, we calculate the Hopf bifurcation threshold th and the pure imaginary eigenvalue ImA corresponding to an 



oscillatory profile instability. For £> = 1, and for a few values of A, in Fig. 19(a) we plot th vs. r associated with 
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Figure 18. For e = 0.02, k — A, D = 1, we plot the phase diagram A vs. r for a four-spot quasi- equilibrium pattern 
in a square with spots centered along the diagonals of the unit square at a distance r from the midpoint of the square. 



the in-phase oscillation eigenvector (1,1,1,1)-^ (heavy sohd curves) and with the out-ol-phase oscillation eigenvector 
(1,-1,1,-1) ( thin solid curves). The dotted portions of the heavy solid curves correspond to where a competition 
instability occurs. From Fig. 19(a)[ we observe that th is set by the synchronous oscillatory instability and that, at 
each fixed r, th increases with A. In Fig. |19(b)| we plot ImA vs. r at the Hopf threshold. From this figure, we observe 
that as r decreases, the pure imaginary eigenvalue for the out-of-phase oscillation eigenvector (1,-1,1,-1)"^ (thin 
solid curve) approaches zero, which implies the onset of a competition instability with eigenvalue A = 0. 




0.1 0.2 0.3 0.4 0.5 0.6 
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(a) Th ^s.. r 
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(b) ImA vs. r 



Figure 19. Fix e = 0.02, /c = 4, and D — 1. For a few values of A we plot the Hopf bifurcation threshold tr 
(left figure ) and pure imaginary eigenvalue A ( right figure ) for a four-spot pattern in the unit square. The in-phase 
oscillation eigenvector (1, 1, 1, 1)"^ corresponds to the heavy solid curves, while the out-of-phase oscillation eigenvector 
(1,-1,1,-1) corresponds to the thin solid curves. From top to bottom, the three heavy solid curves at r = 0.65 
correspond to A — 1.0, 0.9, 0.8. 



We now perform a few numerical experiments to test the predictions of the asymptotic theory. In Experiments 
6.1-6.3 we fix 13 = 1 and e = 0.02, and we vary A, r, and the initial distance r(0) of the spots from the midpoint of 
the square. In these experiments we will focus on competition and oscillatory instabilities of the four-spot pattern. 

Experiment 6.1: (Slow Drift of the Spots): We fix r(0) = 3/(5a/2) and A = 1, corresponding to Regime C of the 
phase diagram in Fig. [TH] We set t — 1, which is below the numerically computed Hopf bifurcation threshold of 
TH = 44.8 at r(0). Thus, the initial four-spot pattern is predicted to be stable to all three instabilities. From the 
dynamics in Principal Result 3.1 we can verify that the four spots drift slowly towards the midpoint of the square. 
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and reach an equilibrium state with r « 1/(2 a/2) for t ^ 1. For t = 1, there is no dynamically-triggered oscillatory 
profile instability. The full numerical results in Fig. 1201 confirm this prediction of the asymptotic theory. 



^mm^m ''■^■H ''■^■H 
o.5^^^H o.5^^^H 0-5|^H O-^I^H ^-^I^H 

0*^™* 0^^^^ 0^^^^ 0^^^^ 0^^^™ 

0.5 1 0.5 1 0.5 1 0.5 1 0.5 1 

t=1 t=101 t = 251 t = 451 t = 701 



Figure 20. Experiment 6.1: Fix e = 0.02, k = A, D = 1, A = 1.0, and t = 1. Starting from r(0) = 3/(5\/2), we 
observe the slow drift of the four spots towards their equilibrium locations within the unit square. At t — 701, they 
are at (approximately) (0.75,0.75), (0.25,0.75), (0.25,0.25) and (0.75,0.25), namely r = 0.5/V2. 

Experiment 6.2: ( A Competition Instability) : Let r = l/(2\/2), corresponding (roughly) to the equilibrium distance 
from the midpoint for the four-spot pattern. We choose A = 0.6 corresponding to Regime /3 in the phase diagram 
of Fig. [iSl In addition, we take r = 1 for which t < th- Therefore, we predict that the initial four-spot pattern is 
unstable to a competition instability This is confirmed by the full numerical results computed from the GS model 
(|1.1[) shown in Fig. [5T] In this figure, we observe that the two spots at (0.25, 0.75) and (0.75, 0.25) decay very fast, 
and that after t — 20, these two spots are annihilated. The spot amplitudes as a function of t are shown in Fig. [22l 
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Figure 21. Experiment 6.2: Fix e — 0.02, D = 1, A = 0.6, and t = 1. This plot shows a competition instability for 
the four-spot pattern with spots initially located at a distance r = l/(2\/2) from the midpoint of the square. 

Experiment 6.3: (An Oscillatory Instability): We fix r = l/(2-\/2) at t = and choose A — 0.8. This choice corre- 
sponds to Regime C, in Fig. 1181 With these parameter values, the numerically computed Hopf bifurcation threshold 
predicted by the asymptotic theory is t^t = 11.0 when f = 0. To confirm this prediction, we compute full numerical 
solutions to the GS model (11.11) for r = 10, 11, 12. For these three choices of t, in Fig. [23] we plot the numerically 
computed amplitude Vm — v„i{t) of the spot at (0.75,0.75). Our full numerical computations show that all the spot 
amplitudes coincide with the one in Fig. [23l This synchronous oscillatory instability was predicted by the asymptotic 
theory. Our numerical results show that an unstable oscillation develops on the range 11 < r < 12, which was 
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Figure 22. Experiment 6.2: Fix e = 0.02, r(0) = l/(2^/2), D^l, 0.6, and r = 1. We plot the amplitude of the 
two spots at (0.25,0.25) (heavy solid curve) and (0.25,0.75) (dashed curve). 

closely predicted by our asymptotic theory. Our numerical results show that this unstable oscillation leads to the 
annihilation of the spot (not shown). As a result, we conjecture that the Hopf bifurcation is subcritical. 
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Figure 23. Experiment 6.3: Fix e = 0.02, r(0) = l/{2^/2), D — 1, and A — 0.8. For the initial four-spot pattern in 
a square, we plot the numerically computed spot amplitude = Vm{t) for the spot at (0.75, 0.75) for t = 10, 11, 12. 
When T = 12 > th , the oscillation is unstable. 



6.2 A One-Spot Pattern in the Unit Disk 

Let VL be the unit disk, and consider a one-spot quasi-equilibrium solution centered at xi G fi, with r = |xi|. From 
()2.22|) and (|2.5|) . the source strength S for this spot is determined in terms of A by 

S + 2t:uRi,iS + vx{S) , A^AeHv^/D), v = -l/\ne. (6.3) 

Here the regular part = i?(xi;xi) of G, which depends on r = |xi|, is calculated from (jA.Sp of Appendix A. 
For this case, the spot dynamics from Principal Result 3.1 shows that the spot will drift with speed 0{e'^) along a 
ray towards the center of the disk. 

For a one-spot solution centered at the origin, the spot self-replication threshold A — As {D) is given by 

As{D) = [S2 (1 -I- 27ri.i?i,i) + (S2)] , (6.4) 

where ^(^2) ~ —1.79 and E2 ~ 4.31. Here = i?i,i(0;0) is given explicitly in (j2.15p . Moreover, the existence 
threshold for a one-spot solution centered at the origin is obtained by replacing E2 in (16. 4p with S, and then computing 
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the minimum value, denoted hy Af = Af{D), of the resulting expression with respect to S. The self-replication and 
existence thresholds for a one-spot solution centered at the origin are shown in the phase diagram of Fig. 1241 
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Figure 24. Plot of the phase diagram A vs. D for e — 0.02 for a one-spot solution centered at the origin in the unit 
disk. It is only in Regime a that spot self-replication can occur. The solid curve is the spot self-replication threshold 
As{D) of (|6.4p . while the lower dotted curve is the existence threshold Af{D). In Regime C a one-spot solution does 
not exist. In Regimes j3 and 7 an oscillatory instability for a spot at the origin occurs only if t > th{0). In Regime 
7, '''^(0) > 0, while in Regime 13, t'^{Q) < 0, where th — TH(r). Thus, only in Regime 7 can we find a value of t 
and an initial spot location for which a dynamical oscillatory instability is triggered for a spot that slowly drifts to 
the center of the unit disk. 
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Figure 25. A One-spot solution in the unit disk at a distance r — |xi| from the origin. Plots of th vs. r for different A 
and D. (a) Fix A — 0.16. The solid curves are for D — A, 5, arranged from lower to upper y— intercepts, respectively. 
The dashed curves are for D — 6,7,8, arranged from upper to lower y-intercepts, respectively, (b) Fix A = 0.18. 
The solid curves are for D ~ 2,3 arranged from lower to upper y— intercepts, respectively. The dashed curves are for 
D — 4.3,5,6, arranged from upper to lower y— intercepts, respectively. 

Next, we study oscillatory instabilities of the spot profile. For each r — |xi| in < r < 1, we suppose that 
A > Af{D,r), where Af — Af{D,r) is the existence threshold for a one-spot quasi-equilibrium solution centered at 



Xi. In Fig. 25(a) we plot the numerically computed Hopf bifurcation threshold th versus r for fixed A — 0.16, but 



for different values of D. The computations were done using the global eigenvalue problem of Principal Result 4.3. 



For D — 4 and D — 5 we observe from Fig. 25(a) that the maximum of th = Tnir) is obtained at the equilibrium 



location r = 0, and that th decreases monotonically as r increases. In contrast, for D — 6,7,8 Fig. 25(a) shows 
that the curves of th = Tnir) are convex near r = 0, so that th has a local minimum at r = 0. Therefore, when 
D is sufficiently large, we conclude that we can obtain a dynamically-triggered oscillatory instability in the spot 
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amplitude. To illustrate this, suppose that D — 8. Then, we calculate th{0) « 3.73 and r/f (0.612) « 3.88. Suppose 
that we take r — 3.8 with the initial spot location at time t — at r — 0.612. Then, since r < r//(0.612), the 
spot is stable at i = 0. However, since the motion of the spot is towards the origin and t > th{0), it follows that 
a dynamically-triggered oscillatory profile instability will occur before the spot reaches the center of the disk. A 



qualitatively similar scenario occurs for other values of A. In particular, for A = 0.18, in Fig. 25(b) we plot th versus 



r for various fixed values of D. For the values D — 4.3,5,6, we observe that Tn^r) has a local minimum at the 
equilibrium location r = 0. Thus, for these values of D, we can choose a value of r and an initial spot location that 
will lead to a dynamically-triggered oscillatory instability. 

For a spot centered at the origin of the unit disk, and with e = 0.02, in Fig. [26] we plot the Hopf bifurcation 
threshold th vs. D ioi A ^ 0.16, 0.21. We observe that each curve is not monotone in D, and has a local 
maximum at some value of D, with th decreasing as D increases. 
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Figure 26. One-spot solution centered at the origin in the unit disk. The Hop} bifurcation threshold Th vs. D is plotted 
for A = 0.16, . . . ,0.21. This shows that th{D) has a maximum at some critical value of D. 



In Fig. [Ml we plot a phase diagram in the A versus D parameter plane for e ~ 0.02 for a one-spot solution centered 
at the origin in the unit disk. From this figure, there are four distinct regions with different solution behavior. In 
Regime a the spot at the origin will undergo self-replication. In Regime C, there is no one-spot solution centered at the 
origin. In Regime /3 the one-spot solution undergoes a Hopf bifurcation if r is large enough and the Hopf bifurcation 
threshold th — Tnir), shown in Fig.[25l satisfies t-^(O) < 0, and so has a local maximum at r = 0. Finally, in Regime 
7, the one-spot solution has a Hopf bifurcation if r is large enough, but now the Hopf bifurcation threshold satisfies 
t'/j{0) > 0, and so has a local minimum at r = 0. Therefore, in Regime 7 one can choose an initial spot location inside 
the unit disk, together with a value of r, for which a dynamic oscillatory instability of a one-spot quasi-equilibrium 
solution will be triggered as the spot drifts towards its equilibrium location at the center of the unit disk. 

We now argue that the solid curve A — As{D) in Fig.[24]not only corresponds to the threshold condition for the self- 
replication of a spot centered at the origin, but also corresponds to the threshold for a dynamically-triggered spot self- 
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replication event for a one-spot pattern in the unit disk starting from the initial spot location at radius r = |xi | with 
< r < 1. To show this, we differentiate (j6.3|) with respect to r to obtain that dS/dr = —2tivS {dRi^i/dr) + 0{v'^). 
It is readily verified numerically by using (lA.Sp of Appendix A that Ri^i has a global minimum at r = with 
an increasing function of r on < r < 1, and Ri,i — > +oo as r — > Therefore, we have dS/dr < on < r < 1, 
so that S' is a monotone decreasing function of r for any fixed values of A and D. Moreover, by difi^erentiating (I6.3p . 
we readily observe that S* is a monotone increasing function of A at each fixed r and D. Therefore, it follows that 
the threshold curve for a dynamically-triggered spot self-replication event is obtained by setting S — a-nd r — 
in (|6.3p to obtain As{D) as given in (j6.4p . Therefore, in Regime a of Fig. [24] it follows that we can choose an initial 
spot location and a value of A for which a dynamic spot self-replication instability will occur for a spot which slowly 
drifts towards the center of the unit disk. 

Experiment 6.4-' (One-spot solution in the unit disk: Comparison with full numerical simulations): For a one-spot 
solution centered at the origin of the unit disk we now compare our numerically computed Hopf bifurcation threshold, 
as obtained from our global eigenvalue problem, with full numerical results as computed from the full GS model (|l.ip . 
For three values of r, in the top row of Fig. [571 we plot Vm = v{0, t) versus t for the GS parameters e = 0.02, A = 0.16, 
and D — 4. This figure shows that an unstable oscillation occurs near r — 3.4, which agrees rather closely with the 
Hopf bifurcation threshold th ~ 3.5 computed from our global eigenvalue formulation of Principal Result 4.3 (see 
Fig. 25(a) I. We remark that if we changed the parameters to A — 0.18 and D = 6, then th ~ 10.3 (see Fig. |25(b)] ), 
as predicted by our global eigenvalue problem. This compares rather closely with the result r w 10 computed from 
the full GS model (II. ip observed in the bottom row of Fig. [57] 




(d) T = 9.0 (e) r = 9.5 (f) r = 10.0 



Figure 27. Experiment 6.4-' One-spot solution with spot at the center of the unit disk. We fix e — 0.02. Top Row: For 
A = 0.16 and D = A, the spot amplitude Vm vs. t is plotted for (a) t — 3.2; (h) r = 3.3; (c) t — 3.4. Bottom Row: 
For A — 0.18 and D — 6, the spot amplitude Vm vs. t is plotted for (d) t — 9.0; (e) t — 9.5; (f) t = 10.0. 
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Table 1. Equilibrium ring radius r^ for a two-spot pattern of the form (j6.5p in the unit disk for different values of D. 

6.3 A One-Ring Pattern of X-Spots in the Unit Disk 

Next, we consider the special case when k spots are equally-spaced on a ring of radius r centered at the origin of the 
unit disk. Representing points as complex numbers in the unit disk, the centers of the spots are at 



J = 1, . 



1. 



(6.5) 



For this symmetric arrangement of spots, the matrix Q in (12.101) is circulant symmetric and is a function of the ring 
radius r. Hence, — 6e with e = (1, . . . , 1)"^, where 9 — pk{r)/k with Pk{r) = J2i=i ^j=i ^iJ- I? ^ 1, we 
can analytically determine a two-term expansion for Pk{r) by using (|2.12[) . the simple explicit form (|A.9|) for the 
Neumann Green's function, and Proposition 4.3 of |32j . In this way, for I? ^ 1, we get 

k k 

kr^-^)+\n{l - r- 



pu{r)^yy^g.^ 



i=i j=i 



k 

"2^ 



\n{kr''-^)+\n(l-r^'')-r^k 



1=1 j=i 

Since Q is circulant, (|2.10l) has a solution for which the source strengths Sj for j = 1, . 
Sc- From p.22p of ^ it follows that Sc satisfies the scalar nonlinear algebraic equation 



3fc 

T 



0{D-^) . (6.6) 
, k have a common value 

k ) v^D me 

To determine the dynamics of the k spots, we calculate in p.Sp as f,- — irk^ -^Scp'kir) e^^J/fc. Then, the dynamics 
p.7p of Principal Result 3.1 shows that all the spots remain on a ring of radius r{t), which satisfies the ODE 



A = S4l + ^Pk{r)]+,^x{Sc), 



A = 



A, 



-1 



'dt ^^i^c)ScPk(r)- 



(6.8) 



This ODE is coupled to the nonlinear algebraic equation (|6.7I) for Sc in terms of r. 

The equilibrium ring radius r^ for (|6.8p satisfies p'j^{r) = 0, and is independent of A. For a two-spot pattern, 
where k ~ 2, then r^ is given in Table [1] for various values of D. There is no simple explicit formula for Pk{r) when 
D = 0(1). However, for D ^ 1, Pk{r) is given asymptotically in ()6.6p . By differentiating this asymptotic result, it 
readily follows that re, for ^ 1, is the unique root on < re < 1 of [fc — l]/(2fc) — = r^''/(l — r"^^) for any 
fc > 2. It is also readily shown from (|6.6p that p'l{re) > 0, with p^(r) < for r < re and p'i^{r) > for r > r^. Hence, 
for D ^ 1, Pk(r) attains its global minimum value at the equilibrium radius re. Therefore, since j{Sc) > 0, then r^ 
is a stable equilibrium point for the ODE (|6.8p when D ^ 1 and, moreover, r — !■ re for any initial point r(0), with 
< r(0) < 1, as i ^ oo. By differentiating ()6.7p for the source strength, this condition on Pk{r) also implies that 
Sc{r) has a maximum value at r — r^ when I? 3> 1. 

Spot Self-Replication Instabilities: We first discuss spot self-replication instabilities for equally spaced spots on 
a ring. For e = 0.02 and D — 3.912, in Fig. [28] we plot the minimum value Af = Af{r) of A for which a quasi- 
equilibrium ring pattern exists at each fixed r for /c = 3, 4, 5 equally spaced spots on a ring. These plots are obtained 
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by determining the minimum point of the A versus Sc curve in (|6.7p . In these plots we also show the spot-splitting 
threshold A = A(r) obtained by setting Sc = ^2~ 4.31 and xi^^) ~ -1.79 in (|6Jl) . 
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Figure 28. For e = 0.02 and D = 3.912 we plot the existence threshold Af (heavy solid curves) and the spot- splitting 
threshold As (solid curves) as a function of the ring radius r for either k — 3 (left figure), k — A (middle figure), or 
k = 5 ( right figure ) equally spaced spots on a ring. 



In Fig. |29(a)[ we fix D = 3.912 and plot Sc versus r for fc = 3, = 30, for fc = 4, = 40, and for fc = 5, = 48. For 
these three patterns, the equilibrium states are, respectively, {rei,Scei) — (0.55,4.57), (re2,S'ce2) — (0.60,4.70), and 
(re3,S'ce3) = (0.63,4.58). These points are the circular dots in Fig. 29(a) The square dots in this figure correspond 
to the spot-replicating threshold S2 ~ 4.31, which occurs at rgi — 0.17, rs2 = 0.14, and rgj, = 0.22. Any portion of 
these curves above the square dots correspond to ring radii where simultaneous spot-splitting will occur. 




CO 




(a) 5c 



(b) Sc 



Figure 29. Fix D ~ 3.912. (a) Sc vs. r with fc = 3, .4 = 30 (dotted curve), fc = 4, ^ = 40 (solid curve) and 
fc = 5, .A = 48 (heavy solid curve), (b) Fix k ~ A, we plot Sc vs. r with ^ = 36 (dotted curve), A = 37.5 (solid 
curve), and A ~ AO (heavy solid curve). The square dots indicate .spot-. splitting thresholds 5c = S2 ~ 4.31, and the 
circular dots, where Sc achieves its maxima, denote equilibrium ring radii. 



In Fig. |29(b)[ we fix fc 4 and D = 3.912 and plot Sc versus r A ^ 36, A ^ 37.5, and A = 40. The equilibria 
are at {re4,Scc4) = (0.60,4.21), {rc5,Sce5) = (0.60,4.40), and irc2,Sce2) = (0.60,4.70), respectively. These points 
correspond to maxima of Sc{r). In this figure the spot-splitting thresholds are the square dots. For A = 36 and fc = 4, 
we observe that since Scd — 4.21 < S2 ~ 4.31 when r — rc, then Sc < ^2 for r 7^ r^, and hence this four-spot pattern 
is stable to spot self-replication for all values of the ring radius. In contrast, consider the solid curve in Fig. |29(b)| 
for A ~ 37.5 and fc = 4. Then, if the initial ring radius r(0) is in the interval between the two square dots in this 
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figure, i.e. 0.35 < r(0) < 0.84, we predict that all four spots will begin to split simultaneously starting at t = 0. 
Similarly, from the heavy solid curve for ^ = 40 and k — A in Fig. |29(b)[ we predict that the four spots will split 
simultaneously when r(0) > rs2 = 0.14. This prediction is confirmed in the full numerical results below. 

The plots of As — As{r) and Sc — Sc{r) in Fig. [28] and Fig. [29l respectively, also clearly show the possibility of a 
dynamically-triggered spot-splitting instability. For instance, consider the solid curve in Fig. [29(b)] for A = 37.5 and 
k — 4: and an initial ring radius r(0) = 0.2. This initial point is below the spot-splitting threshold shown in Fig. 28(b)] 
and so the initial pattern is stable to spot-splitting. However, eventually as t increases the ring radius will cross the 
threshold value r sa 0.35 for spot-splitting, and a dynamically-triggered simultaneous spot-splitting event will occur. 

This possibility of a dynamically-triggered spot-splitting instability for equally spaced spots on a ring for the GS 
model is in distinct contrast to the behavior found in '38' for the Schnakenburg model. For this related RD model, 
the common source strength for a pattern of k equally-spaced spots on a ring is independent of the ring radius r. 
This precludes the existence of a dynamically-triggered spot-splitting instability for the Schakenburg model. 

Experiment 6.5: (One-ring pattern and spot- splitting when D 3> 0{1)): For = 40 and D — 3.912, we consider an 
initial four-spot pattern with equally-spaced spots on an initial ring of radius r = 0.5 > rs2 at time t = Q. Since 
Sc ~ 4.69 at i = 0, as computed from (|6.7|) . our asymptotic theory predicts that all four spots will begin to split 
simultaneously a.t t — 0. The full numerical results for v computed from the GS model (|l.ll) . as shown in Fig. 1501 
confirm this prediction, and also show that the direction of splitting is perpendicular to the direction of spot motion, 
as predicted by Principal Result 3.2. The splitting process generates an equilibrium pattern of eight equally-spaced 
spots on a ring. 




Figure 30. Experiment 6.5: Fix D — 3.912, A = 40, t — \, and e — 0.02. Consider an initial four-spot pattern with 
equally- spaced spots on a ring of radius r = 0.5. The initial common source strength is Sc ~ 4.69. From left to right 
we plot V at times t — 4.6, 70, 93, 381. All spots undergo .splitting, and the dynamics leads to an equilibrium eight-spot 
pattern on a ring. 



Phase Diagrams for Existence, Self- Replication, and Competition: Next, we use our global eigenvalue prob- 



lem of Principal Result 4.3 of i l4.2l to compute the competition instability threshold as a function of the ring radius 
r for various values of the parameters. Our results are shown in terms of phase diagrams of A versus r. 

For e — 0.02, in Fig. [31] we plot the phase diagram A versus r for D = 0.2, D = 1.0, and D = 5.0, for a two-spot 
quasi-equilibrium solution, showing the thresholds for the existence of a quasi-equilibrium pattern, for a competition 
instability, and for a spot self-replication instability. For each fixed r, the existence thresholds in these figures were 
computed by determining the minimum of the curve A versus Sc in (|6.7I) , while the spot self- replication thresholds 
were obtained by substituting S'c = S2 ~ 4.31 and ~ —1.79 in ()6.7p . In addition, the competition instability 

thresholds were computed from (|6.9p . as we discuss below. 

The solution behavior in the four distinct parameter regimes of Fig. ]3T] is described in the caption of Fig. ]3I] As 
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expected, the subfigure for D = 0.2 on the range < r < 0.5 is quahtatively similar to the phase diagram shown in 



Fig. 14(a) for a two-spot solution on the infinite plane. In Fig. |32]we show similar phase diagrams for D = 0.2 and 
e = 0.02, but for fc = 4, fc = 8, and k = 16, equally-spaced spots on a ring of radius r. 
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Figure 31. The phase diagram A vs. r for a two-spot quasi- equilibrium solution on a ring of radius r in the unit disk 
when e = 0.02. The solid curve is the existence threshold Af, the dotted curve is the competition instability threshold; 
the heavy solid curve is the spot self -replication threshold A^. Regime a: no two-spot solution; Regime /3: the two-spot 
solution is unstable to competition; Regime ^ the two-spot solution is unstable to a oscillation if t > T/f(r); Regime 
9: the two-spot solution is unstable to spot self-replication. Subfigures from left to right:: (a) D — 0.2; (b) D — 1; (c) 
D = 5. 
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Figure 32. The phase diagram A vs. r for a k-spot quasi- equilibrium solution on a ring of radius r in the unit disk 
with e =^ 0.02 and D = 0.2. The thin solid curves are the existence threshold Af for the quasi- equilibrium solution, 
the dotted curves are the critical values of A for a competition instability; the heavy solid curves are the spot self- 
replication threshold Ag. In Regime a the quasi- equilibrium solution does not exist. In Regime (3 the quasi- equilibrium 
solution exists but is unstable to competition. In Regime ( the solution is unstable to an oscillatory instability if t 
exceeds a Hopf bifurcation threshold th{t). In Regime 9 the solution is unstable to spot self- replication. Subfigures 
from left to right: (a) k = A; (b) k = 8; (c) k = 16. 



From §4.2l and the result ()4.20p for the spectrum of the A-dependent circulant Green's matrix Q\, it follows for /c = 4 
that the eigenvectors of Qx arc vi ~ (1,1,1,1)-^, V2 = (1,0,-1,0)^, V3 = (1,-1,1,-1)-^, and V4 = (0,1,0,-1)"^. 
Note that the last three of these vectors have components with different signs. Our computational results show that 
the competition instability threshold is set by the sign- fluctuating instability V3, which has a larger competition 
instability regime in the phase diagram of A versus r than does either V2 or V4. There is no competition instability 
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threshold associated with vi. Similarly for fc = 8 and fc = 16, the instability associated with the eigenvector of the 
form V = (1, —1, 1, —1, . . . , 1, —1)"^ is always the dominant competition mechanism. In contrast, as we discuss below, 
our numerical results show that the Hopf bifurcation threshold th for an oscillatory profile instability is set by the 
eigenvector Vi , which corresponds to a synchronous oscillatory instability. 

Since the sign-fluctuating eigenvector determines the competition instability threshold in Fig. [3T] and Fig. |32l we 
get from (|4.24p that, for fc even, this threshold is obtained by numerically solving the simple coupled system 
fc 

X'(5c) + 27ry A^^^[Sa{l + 27riy0) + i^xiSc)] , ^ ^ , (6-9 a) 

WD In e 



m— 1 

fc 



a,„ , fli = Rii ; Qj =Gij , j = 2,...,k. (6.9 b) 



m— 1 



Here 6 is the eigenvalue of G with eigenvector e = (1, . . . , 1)^. Numerical values for Ru and Gij, for j = 2, . . . , fc, 
are obtained from the reduced-wave Green's function for the unit disk given in (|A.8p of Appendix A. 

Oscillatory and Competition Instabilities: Next, we use our global eigenvalue problem to compute Hopf bifur- 
cation thresholds th for an oscillatory instability of fc equally-spaced spots on a ring. Such an instability is the only 
instability mechanism in each Regime C of Fig. |32]and Fig. |3T] Since this instability is not readily depicted in a phase 
diagram of the type shown in Fig. 1321 and Fig. 1311 we now discuss it some detail. 

To focus the discussion, we will only consider an equally-spaced two-spot pattern on a ring of radius r in the 
unit disk. We show that both dynamically-triggered oscillatory and competition instabilities are possible for different 
ranges of the GS parameters, and we study the spectrum of the linearization in some detail. 

For A = 0.26, fc = 2, and e = 0.02, the Hopf bifurcation threshold th versus the ring radius r, corresponding to 
Vi — (1,1)^, is shown by the solid curves in Fig. [33l The dotted segments on these curves is where a competition 
instability occurs. In Fig. [33(a)] we plot TH{r) for D = 0.8, 1.0, 1.3, 1.5, while in Fig. [33(b)] we plot rff(r) for D = 
2, 3, 4, 5. From the convexity of th versus r in Fig. [33(b)[ we conclude that a dynamically-triggered oscillatory profile 
instability can occur for a two-spot pattern when A — 0.26 and D > 2. We also observe that when _D = 5, a 
competition instability occurs when r 0.568. Therefore, since the equilibrium ring radius is r^ ~ 0.455 for D = 5 
(see Table [ij, we conclude that there can be a dynamically-triggered competition instability for this value of D. To 
illustrate this, let D = 5 and r = 1 and suppose that we have an initial two-spot pattern with initial value r — 0.7 at 
t = 0. Then, as the two spots move slowly towards each other, a dynamically-triggered competition instability will 
be initiated before the ring reaches its equilibrium radius at r^ ~ 0.455. 

Next we fix A = 0.26 and D ~ 3. In Fig. 34(a) we plot the Hopf bifurcation thresholds thi and th2 for in-phase 
and out-of-phase spot amplitude oscillations, respectively. Since thi < th2, the oscillatory instability triggered as r 
increases is a synchronous spot amplitude oscillation. The dotted portion on thi, representing the continuation of 
the heavy solid curve in Fig. 34(a) shows where a competition instability occurs for any r > 0. In Fig. [34(b)[ we 
fix two spots on an equilibrium ring of radius r^ = 0.45483, and we plot the path of the complex conjugate pair of 
eigenvalues in the complex plane for the range r G [3.7, 36]. We note that as r increases, the real parts of both pairs 
of eigenvalues increase. The eigenvalue for in-phase oscillations Vi (1, 1)-^, which is given by the solid curve in 
Fig. [34(b)] enters the right half-plane at r w 11.5. Note that the imaginary part of the eigenvalue for the out-of-phase 
oscillation V2 = (1,-1)"^ becomes very close to the negative real axis when r < 3.7. 

Next, for e — 0.02 and A = 0.26, we treat D as the bifurcation parameter and we relate our results with those 
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Figure 33. Two equally- spaced spots on a ring of radius r in the unit disk for e ~ 0.02; (a) For A = 0.26, the Hopf 
bifurcation threshold th vs. r is plotted for D = 0.8, 1.0, 1.3, 1.5, with the lower curves corresponding to smaller 
values of D. The solid curves correspond to regions in r where an oscillatory profile instability occurs, and the dotted 
portions of these curves correspond to where a competition instability occurs for any r > 0. (b) th vs. r for the larger 
values of D given by D = 2,3,4,5. The lower curves at r = 0.6 correspond to larger values of D. 
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Figure 34. Two equally- spaced spots on a ring of radius r in the unit disk. Fix A = 0.26, e ~ 0.02, and Z? = 3. 
(a) Th vs. r for in-phase spot amplitude oscillation vi = (1, 1)-^ (lower heavy solid curve), and th2 for the out-of- 
phase oscillation V2 = (1,-1)"^ (upper solid curve). The dotted portions on thi are where a competition instability 
occurs for any r > 0. (b) For the equilibrium ring radius r^ — 0.45483, we plot A in the complex plane on the range 
T G [3.7,36]. The solid curve isforvi = (1,1)-^, while the dotted curve is for V2 = (1,-1)^. As t increases (direction 
of the arrows), i?e(A) increases. 



obtained from the NLEP theory of |62) . In Theorem 2.3 of [62], as summarized in Appendix B, the leading-order-in-r/ 
NLEP analysis proves that if £) = and Lq < (2i]o+k)'^ ' ^^^^ *^he small solution u^,v^ of (|l.ll) is stable for any 

T sufficiently small. In addition, the NLEP theory proves that if Lq > -(2ri^+k)^' then this solution is unstable for any 
T > 0. From (jB.7|) . we recall that Lo ~ hm£_j.o j^j" and rjo = lini£_j.o '^here bo — w^pdp sa 4.9347 and 

w{p) satisfies (|B.1[) . With = tt and h' = — 1/lne, and for the other GS parameter values as given, this stability 
bound is 1.158 < D < 3.303. For this range of D, the leading-order-in-j/ NLEP theory predicts that this solution 
is stable for t sufficiently small. Outside of this bound for D, the NLEP theory of |62j predicts that a competition 
instability occurs for any r > 0. This conclusion from NLEP theory is independent of the ring radius r. 



For e = 0.02 and A = 0.26, in Fig. 35(a) we plot Dc versus r by the solid curve, where Dc denotes the critical 
value of D where a competition instability is initiated, as computed from (16.91) . To the left of this curve, oriented 
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Figure 35. Two equally- spaced spots on a ring of radius r in the unit disk. Fix A — 0.26 and e = 0.02. (a) Plot of 
Dc vs. r (solid curve), where Dc is the critical value of D at which a competition instability is initiated. The region 
1.1583 < D < 3.3030 between the two horizontal dotted lines is where the leading- order-in-v NLEP theory predicts 
that no competition instability occurs, (b) For r = 30 and r ~ 0.3, a plot of the spectrum in the complex plane is 
shown when D varies on the range [0.8, 3.1], which corresponds to taking a vertical slice in Fig. 35(a) at r = 0.3. The 
solid path in the spectrum corresponds to V2 = (1, —1)"^, while the the dotted path is for vi = (1, 1)-'". The arrows 
show the path as D increases. 
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Figure 36. Two equally- spaced spots on a ring of radius r — 0.3 in the unit disk. Fix A — 0.26, e — 0.02, and r — 30, 
and consider the sign-fluctuating eigenvector V2 = (1, — 1)"^- We plot the eigenvalues of the global eigenvalue problem 
for two ranges of D where they are real. Left figure: D G [0.8,0.845]. Right figure: D G [2.95,3.10]. 



with respect to the direction of increasing £>, we predict that a competition instabihty occurs for any r > 0, while 
to the right of this curve we predict that the solution is stable when r is small enough. The two horizontal dotted 
lines with Di = 1.158 and D2 = 3.303 in this figure bound a region of stability as predicted by the leading-order-in-j/ 
NLEP theory of [62] . This prediction from NLEP theory does not capture the dependence of the stability threshold 
on the ring radius r. Our stability formulation, which accounts for all orders in = — 1/ Ine, shows that Dc = Dc{r). 
For £ = 0.02, v is certainly not very small, and this dependence of the threshold on r is very significant. 



In Fig. 35(b) for A = 0.26, r = 30, and r = 0.3, we plot the spectrum in the complex plane as D varies over the 
range [0.8,3.0]. This range of D corresponds to taking a vertical slice at r = 0.3 in Fig. 35(a) that cuts across the 
stability boundary Dc{r) at two values of D. The arrows in Fig. 35(b) indicate the direction of the path of eigenvalues 
as D is increased from D = 0.8. The loop of spectra, depicted by the solid curves in Fig. |35(b)| shows the eigenvalues 
associated with the eigenvector V2 ~ (1,-1)"^. From Fig. 35(a) we predict instability when D is near D = 0.8 or 
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when D is near D = 3.0. This behavior is shown by the closed loop of spectra in the complex A plane of Fig. |35(b)[ 
and in Fig. |36] where the eigenvalues are plotted for the range of D where they are real. 

More specifically, with regards to the loop of spectra (solid curves) in Fig. 35(b) associated with the V2 = (1, —1)"^ 
eigenvector, our results show that for D = 0.8 there are two real eigenvalues at A w 0.218 and A ~ —0.0224 (the 
negative eigenvalue is not shown in Fig. [35(b)] ). From Fig. [36(a)] we observe that as D is increased above D — 0.8 the 
positive real eigenvalues move along the horizontal axis in different directions, and collide aX D k, 0.845, producing 
a complex conjugate pair of unstable eigenvalues at this value of D. Then, for D w 1.0, the complex conjugate pair 
of eigenvalues enters the stable left half-plane Re(A) < 0. For D = 2.89 the eigenvalues merge onto the negative real 
axis, as shown in Fig. [36(b)] One eigenvalue then enters the right half-plane at _D w 2.95, while the other slides along 
the negative real axis. Therefore, for D > 2.95, there is a positive real eigenvalue associated with the eigenvector 
V2 = (1, —1)"^, which generates a competition instability. 

Alternatively, the dotted curves in Fig. 35(b) correspond to the path of eigenvalues associated with the vi = (1, 1)^ 
eigenvector, representing synchronous oscillatory instabilities. This path depends sensitively on the value chosen for 
T. For T = 30, our computational results show that as D is increased above D = 0.8, the real part of these eigenvalues 
first starts to decrease and then crosses into the stable left half-plane Re(A) < at I? w 1.33. The real part of these 
eigenvalues starts increasing when D w 1.59 and the path then re-enters the unstable right half-plane a.t D — 1.89 
where stability is lost at a Hopf bifurcation corresponding to synchronous oscillations of the two spot amplitudes. 
This path of eigenvalues remain in the unstable right half-plane Re(A) > for D > 1.89. 

For T — 30, A — 0.26, and e = 0.02, we conclude that if 1.33 < D < 1.89, then the two-spot pattern is stable to 
both competition and synchronous oscillatory instabilities. When D > 2.95 there is an unstable real eigenvalue in the 
right half-plane associated with a competition instability as well as a synchronous oscillatory instability associated 
with the complex conjugate pair of eigenvalues. On the range 1.89 < D < 2.95 there is only a synchronous oscillatory 
instability. We remark that if r is taken to be much smaller than our chosen value r = 30, then for D > 2.95 there 
would be only a competition instability from a positive real eigenvalue, and no oscillatory instability would occur. 

Experiment 6.6: (One-ring pattern with two equally- spaced spots on a ring: Simultaneously occurring instabilities): 
For an equally-spaced two-spot pattern on a ring of radius r = 0.3, we now validate the theoretical prediction, based 
on Fig. [35l of having both a competition and oscillatory instability occurring simultaneously for the GS parameter 
values A — 0.26, D — 3, t = 30, and e — 0.02. For these values. Fig. [35(b)[ shows that there is an unstable real 
eigenvalue together with an unstable complex conjugate pair of eigenvalues. In Fig. [37] we plot the maxima of v 
versus t of the two spots, as computed from a full numerical solution of (jl.ip . This figure shows that one of the two 
spots is annihilated before t — 13, consistent with a competition instability. In contrast, the amplitude of the other 
spot begins to oscillate, but the oscillation ceases after the first spot is annihilated. This occurs since after one spot 
has been destroyed, the value r = 30 is below the Hopf bifurcation threshold for a one-spot solution in the unit disk. 



7 Some Asymmetric Spot Patterns in the Unit Square and Disk 

In this section we study spot self-replication instabilities for some spot patterns for which the Green's matrix is 
not circulant. The GS parameters are chosen so that competition and oscillatory instabilities do not occur. The 
asymptotic results for spot self-replication and spot dynamics are favorably compared with full numerical results. 
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Figure 37. Experiment 6.6: Two equally- spaced spots on a ring of radius r = 0.3 in the unit disk. Fix e — 0.02,74 = 
0.26, D — Z, and t = 30. The dot-dashed curve plots the amplitude v„ii vs. t of one the spots, while the solid curve 
shows the amplitude Vm2 vs. t of the other spot. Both a competition and an oscillatory instability occur for the initial 
pattern. 
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Table 2. Experiment 7.1: The columns of Xj, yj, Sj and Vj{0) correspond to the initial condition. The remaining 
data is with regards to the numerical and asymptotic results, respectively, for the equilibrium state. 



7.1 The Unit Square 

We will consider three examples for the unit square [0, 1] x [0, 1] for e = 0.02 and r = 1.0. 

Experiment 7.1: (A three-spot pattern: Slowly drifting spots): Let ^ = 20 and D — 1, and consider an initial three- 
spot pattern with spots equally-spaced on a ring of radius r = 0.2 centered at the midpoint of the unit square. 
The initial spot coordinates at t = are given in Table [2] For this pattern, the matrix Q in ()2.10p is not circulant. 
However, by evaluating the entries of G by using (|A.lip of Appendix A, we compute that the three spots still have 
an initial common source strength Sc ~ 3.71. Since Sc < S2, we predict that there is no spot self-replication initiated 
at t = 0. From the full numerical results shown in Fig. 1381 '^6 observe that all three spots drift slowly outwards, and 
approach equilibrium locations inside fl when t is large. During this evolution, the source strengths never exceed the 
threshold S2, and so there is no dynamically-triggered spot self-replication instability. 

In Fig. we show a very favorable comparison between the spot trajectories, as obtained from the dynamics 
p.7|) of Principal Result 3.1, and the corresponding numerical results computed from (jl.l|) . Since the initial locations 
of the 1^^ and 2^^"^ spots are exactly symmetric, it follows that G1.3 — 6*2,3, ^1.1 = R2.2, and Gi,2 — G2,i in the 
Green's matrix Q. Thus, from Principal Result 3.1, these two spots move at the same speed and their trajectories 
in the a;— direction exactly overlap, as shown in Fig. 39(a) Moreover, from the spot equilibrium condition (j3.8p . we 
can calculate the equilibrium locations Xje,yje, and source strengths Sje for j = 1, ... ,3. These values are given 
in Table [21 where we observe that our equilibrium result closely predicts the final equilibrium state Xjn , yjn , Sj„ 
computed from full numerical solutions of the GS model (|1.1[) . 
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Figure 38. Experiment 7.1: Fix A ^ 20, D ^ \, t = \, and e = 0.02. Consider a three-spot initial pattern with spots 
equally- spaced on a ring of radius r — 0.2 centered at (0.5, 0.5) at time t = in the unit square. The initial common 
source strength is Sc ~ 3.71 (see Table\^. The numerical solution v is shown at various times. 




(a) Xj vs. t 



t 

(b) vs. t 



Figure 39. Experiment 7.1: Fix A = 20, D = 1, t = 1, and e = 0.02. Comparison of the asymptotic dynamics of a 
3— spot pattern (solid curves) with full numerical results (discrete markers) starting from t = 10. (a) Xj vs. t; (b) yj 
vs. t. 



Experiment 7.2: (A three-spot pattern: Spot- splitting and dynamics after splitting): Next, we fix = 20 and D = 1 
and consider an initial three-spot pattern on a ring of radius r — 0.3 centered at (0.4,0.4) in the unit square. These 
GS parameter values are the same as in Experiment 7.1, except that here the initial spot configuration is different. 
The initial data are given in the left side of Tabled Since 5*3 > S2, the asymptotic theory predicts that the 3'"'^ spot 
will undergo splitting beginning aX t — Q. In Fig. 1401 we plot the numerical solution v at different instants in time. 
From this figure we observe that the 3'^'^ spot (the brightest one at t = 1) deforms into a peanut shape ai t = 31, 
and then splits into two spots at t = 46 in a direction perpendicular to its motion. Subsequently, the four spots drift 
slowly to a symmetric four-spot equilibrium state. 

After the splitting event, we use the new 4— spot pattern Xjt, yjt for j = 1, . . . , 4 at t = 61, as given in Table [Sj 
as the initial conditions for the asymptotic DAE system of Principal Result 3.1. At t — 61, the source strengths 
Su, ■ . ■ ,S4t, as given in Table [3l are all below the spot-splitting threshold. In Fig. IHl we show a very favorable 
comparison between the asymptotic trajectories of the spots for t > 61 and the full numerical results computed from 
(jl.l|) . Therefore, Principal Result 3.1 closely predicts the motion of a collection of spots after a spot self-replication 
event. We remark that in Fig. 41(a)[ the x— coordinates of the 1^^ and 2^^*^ spots, as well as the 3'"'^ and 4^^ spots, 
essentially overlap as a result of their almost symmetric locations. 

From the spot equilibrium condition (13. 8p . we calculate the asymptotic result for the symmetric equilibrium 
locations Xje, yje and source strengths Sje, for j = 1, . . . ,4. These results are given in Table |3l and compare almost 
exactly with the results from the full numerical simulations of the GS model (|1.1[) at i = 421 (not shown). 
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j 




yj 


Sj 






yjt 


Sjt 






Sje 


1 


0.25 


0.66 


4.05 


32.54 


0.30 


0.75 


3.44 


0.25 


0.75 


2.87 


2 


0.25 


0.14 


2.37 


35.48 


0.28 


0.20 


3.01 


0.25 


0.25 


2.87 


3 


0.70 


0.40 


4.79 


27.65 


0.76 


0.36 


2.40 


0.75 


0.25 


2.87 


4 










0.77 


0.59 


2.55 


0.75 


0.75 


2.87 



Table 3. Experiment 7.2: The columns of Xj, yj, Sj and Vj{0) correspond to the initial condition. The data Xjt, yjt, 
and Sjt correspond to the initial conditions used for the asymptotic dynamics at t = 61 after the spot- splitting event. 
The final columns are the equilibrium results for the A-spot pattern obtained from the asymptotic theory. 




t=61 t=91 t=121 t=421 

Figure 40. Experiment 7.2: Fix A = 20, D = \, t — \, and e — 0.02. Consider a three-spot initial pattern with spots 
equally- spaced on a ring of radius r = 0.3 centered at (0.4, 0.4) in the unit square. The initial source strengths Sj 
for j = 1 , . . . , 3 are given in Table for which S3 > S2 • The full numerical solution for v is plotted at different 
instants in time. The i^'^ spot undergoes self-replication. Subsequently, all four spots slowly drift towards a symmetric 
equilibrium A— spot pattern. 
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Figure 41. Experiment 7.2: Fix A = 20, D = 1, r = 1, and e — 0.02. Starting from i = 61 after the spot- splitting 
event, and using the initial values as in Table\3^ we compare the dynamics of the 4— spot pattern from the asymptotic 
analysis (solid curves) with the full numerical results (discrete markers), (a) Xj vs. t; (b) yj vs. t. 
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j 




V] 




Vj 


Xjt 


2/jt 


Sjt 






Sje 


1 


0.80 


0.60 


2.82 


36.21 


0.84 


0.58 


3.13 


0.83 


0.61 


2.91 


2 


0.60 


0.80 


2.82 


36.21 


0.58 


0.84 


3.13 


0.61 


0.83 


2.91 


3 


0.40 


0.60 


5.69 


19.07 


0.25 


0.49 


2.93 


0.17 


0.39 


2.91 


4 


0.60 


0.40 


5.69 


19.07 


0.26 


0.66 


2.73 


0.21 


0.79 


3.07 


5 










0.49 


0.25 


2.93 


0.39 


0.17 


2.91 


6 










0.66 


0.26 


2.73 


0.79 


0.21 


3.07 



Table 4. Experiment 7.3: The columns of Xj, yj, Sj and Vj{0) correspond to the initial condition. The data Xjt, Ujt, 
and Sjt correspond to the initial conditions used for the asymptotic dynamics at t = 21 after the two spot- splitting 
events. The final columns are the equilibrium results for the 6-spot pattern. 

Experiment 7.3: (An asymmetric four-spot pattern): We fix = 30 and D = 1, and consider an initial 4— spot 
pattern with equally-spaced spots on a ring of radius 0.2 centered at (0.6, 0.6) in the unit square. The initial spot 
locations and source strengths are given on the left side of Table SI For this case, we find that S3 > S2 and 6*4 > S2, 
so that our asymptotic theory predicts that these two specific spots undergo self-replication events starting at i = 0. 
The full numerical results for this pattern are given in Fig. S^l where we observe that the 3'"'^ and 4^-'^ spots split in a 
direction perpendicular to their motion. The resulting 6— spot pattern, with initial locations Xjt, yjt for j — 1, ... ,6 
at t = 21, as given in Tabled are used as initial conditions for the asymptotic DAE system in Principal Result 3.1. 
These six spots evolve outwards as t increases and eventually approach their equilibrium states Xje , yje , and Sje , for 
j = 1, . . . , 6, as given in Tabled when t ~ 581. The asymptotic result obtained from the spot equilibrium condition 
(|3.8p for the 6— spot equilibrium state is found to compare very favorably with the full numerical results at t — 581. 

A very favorable comparison of analytical (solid curves) and full numerical results (discrete markers) for the spot 
trajectories after i = 21 is shown in Fig.SJ] This agreement further validates our spot equilibrium condition (13. 8p . 




t=45 t=141 t=341 t=581 

Figure 42. Experiment 7.3: Fix A — 30, D — I, t — 1, and e — 0.02. Consider a four-spot initial pattern with spots 
equally spaced on a ring of radius r = 0.2 centered at (0.6,0.6) in the unit square. The initial source strengths Sj for 
j = 1, ... ,4 are given in TableU\ for which S3 > E2 and S4 > T12. The full numerical solution for v is shown at 

— W 

different instants in time. The 3 and 4''"' spots undergo self-replication. Subsequently, all six spots move towards a 
symmetric 6— spot equilibrium pattern. 
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Figure 43. Experiment 7.3: Fix A = 30, D ~ I, t = 1, and e = 0.02. For t > 21, after the two spot- splitting events, 
we compare the asymptotic dynamics of the 6 spots (solid curves) with full numerical results (discrete markers), (a) 
Xj vs. t; (b) yj vs. t. 

7.2 The Unit Disk: A Ring Pattern of Spots with a Center Spot 

In this subsection, we let e = 0.02 and t — 1, and we construct a spot pattern for which fc — 1 spots are equally 
spaced on a ring of radius r, with < r < 1, and with an additional spot located at the center of the unit disk. The 
fc — 1 spots on the ring have a common source strength Sc, while the fc^^ spot at the origin has a source strength 
denoted by Sk- For this pattern, the centers of the spots, written as complex numbers inside the unit disk, are at 

= re^^*^'/(''■"^^ .7 = l,...,fc- 1; = , fc>3. 

We will analyze the dynamics of the spots and the possibility of spot-splitting for this pattern for the limiting case 
-D ^ 1, for which the Green's function G can be approximated by the Neumann Green's function G*^^-'. Recall, that 
G*^^-* has a very simple analytical formula in the unit disk (see Appendix A. 2). As a result, the dynamics of the spots 
and the conditions for spot self-replication can be obtained in a very explicit form. 

For this pattern. Principal Result 3.1 shows that the fc spot at the center of the disk is stationary. In addition, 
by differentiating (|6.6|) . it readily follows that the dynamics of the remaining fc — 1 spots on the ring satisfy 

2 / fc — 1 

Here G(r;0) is the radially symmetric reduced-wave Green's function with singularity at the origin, and pk-i{r) is 
defined in ()6.6|) . with limiting asymptotics for _D ^ 1 also given in (|6.6|) . Since Xj — r e^^^'^/C^^^) for j = 1, . . . , fc — 1, 
then (|7.ip reduces to an ODE for the ring radius given by 

2 J fc- 1 

in terms of Sc — Sc{r) and Sk = Sk{r). From (|2.7|) . we obtain that 5c — Sc{r) and 5*^ = Sk{r) are the solutions to 
the coupled two-component nonlinear algebraic system 
2Triy 



X' ^ -2^£2^(5c) 



,fc- 1 



(7.1) 



dr 
'dt 



-2ne'-f{Sc 



+ SkGr{r;Q) 



(7.2 a) 



A^SAl 



r 



Pk^iir) + lyxiSc) + 2^uG{r- Q)Sk , 



A = 2Tiv{k - l)G(r; 0)5^ + vx{Sk) + Sk{l + 2^uRo^o) , 

(7.2 &) 

where i?o^o is the regular part of G{r; 0) at the origin. The problem (|7.2I) is a DAE system for the ring radius, 
consisting of the ODE ()7.2 a|) coupled to the two nonlinear constraints (j7.2 The core problem for an individual 
spot enters only through the determination of xi^) ^-i^d ^{S). 
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Since pfc_i(r) is only known explicitly when D 3> 1, we use the large D asymptotics (j6.6p . together with 



G(r; Q)^D + G^^\r; 0) = D - ^^Inr + ^ - ^, 



Ro,o 



8tt 



-D > 1. 



(7.3) 



to calculate the equilibrium ring radius. By setting r' = in (17.2 ap . we obtain the following transcendental equation 
for the equilibrium ring radius when Z? 3> 1: 

^2fe-2 



^(1 



(fc-2) 



r\k-l) - (fc- 1) 



1 



0. 



(7.4) 



By solving the three nonlinear algebraic equations in (j7.2 b\i and (|7.4p we obtain equilibrium values rg, See, and Sk 



In Fig. [44(a)] and Fig. |44(b)[ we fix I? = 3.912 and we plot the numerically computed source strengths Sc and Sk, 
respectively, versus the ring radius r for fc = 3,y^ 22, and for fc = 4,^ = 38, as computed from (|7.2 b^ . In this 



computation we used the large D asymptotics of (j6.6p and (j7.3p in (j7.2 In Fig. 44(a) the upper (lower) branch 
of the Sc — Sc{r) curve corresponds to the lower (upper) branch of the Sk — Sk{r) curve shown in Fig. |44(b)[ For 
each pattern, there are two equilibrium ring radii, which are marked by the solid dots. The spot-splitting threshold 
S2 is also marked on the curves in Fig. 44(a) and Fig. |44(b)] by the open dots. We now show how these figures can 
be readily used to predict spot-splitting behavior. 
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Figure 44. Fix D = 3.912, e = 0.02. (a) The solid curves correspond to A — 22. fc = 3 and the heavy solid curves 
correspond to A — 38, fc = 4. The open dots correspond to spot- splitting thresholds, and the solid dots correspond to 
equilibrium ring radii, (b) The heavy solid curves correspond to A = 82, fc = 10 and the solid curves correspond to 
A = 60, fc = 10. Left: Sc vs. r; Right: Sk vs. r. Each upper (lower) branch for Sc{r) corresponds to the lower (upper) 
branch for Sk {r) . 



Localized Spot Patterns in the Two- Dimensional Gray-Scott Model 



55 



Experiment 7.4-' (A one-ring pattern with a center-spot for D ^ 0(1) •' Different types of spot- splitting): We fix 
,4 = 38 and D = 3.912 and consider an initial pattern of one center-spot and tfiree equally-spaced spots on a ring 
of initial radius r = 0.8 at time t = 0. We then predict the dynamical behavior of this initial spot pattern from the 
heavy solid curves in Fig. 44(a) and Fig. |44(b)[ From the data used to plot Fig. |44(b)[ we obtain from the upper 
branch of the Sk versus r curve that Sk — 5.38 > S2 at time t = 0. Correspondingly, from the lower branch of the Sc 



versus r curve in Fig. 44(a) we obtain that 5c < S2 at i = 0. Therefore, we predict that the center-spot will undergo 
spot-splitting starting at t = 0. The other three spots remain on a ring whose radius decreases slowly in time. For 
this parameter set, the full numerical results computed from the GS model (jl.ip . shown in the top row of Fig. I45[ 
confirm the prediction based on the asymptotic theory. 

In contrast, suppose that the initial ring radius is decreased from r = 0.8 to r 0.5. Then, from the heavy 
solid curve in Fig. [44(a)| we now calculate that Sc ~ 4.92 > S2 from the lower branch of the Sc versus r curve. 
Correspondingly, from the upper branch of the Sk versus r curve in Fig. |44(b)| we get S'fe < S2 at t = 0. Thus, for this 
initial ring radius, we predict that the three spots on the ring undergo simultaneous spot-splitting starting at i = 0, 
while the center-spot does not split. The full numerical results shown in the second row of Fig. 251 as computed from 
the GS model (jl.ip . fully support this prediction based on the asymptotic theory. 

We remark that based on the relative locations of the equilibrium ring radius and the spot-splitting threshold in 
Fig. |44(a)| and Fig. |44(b)[ we conclude that it is not possible to obtain a dynamically-triggered spot-self replication 
instability for these parameter values. This more exotic instability is shown in our final example. 




Figure 45. Experiment 7.4-' Fix A = 38, D = 3.912, and e = 0.02. The numerical results for v computed from 
the GS model (jl.ip at different instants in time are shown. The initial condition has three equally- spaced spots 
on a ring together with a center-spot. Top row: initial ring radius r ~ 0.8; Plots of v from left to right at times 
t = 4.6, 46, 102, 298. Bottom row: initial ring radius r = 0.5; Plots of v from left to right at times t — 4.6, 56, 140, 298. 



Experiment 7.5: (A one-ring pattern with a center-spot for D 3> 0{1): A dynamically-triggered instability): In 
Fig. |44(c)| and Fig. |44(d)[ we fix fc = 10 and D = 3.912 and plot Sc and Sk versus r ioi A = 82 and A = 60. In 
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these figures, each sohd dot and empty circle denotes an equilibrium ring radius and the spot-splitting threshold, 
respectively. 

Consider first the solid curve in Fig. 44(c) and Fig. |44(d)| corresponding to .4 = 60. For this value of A it follows 
that only the center-spot can split and the minimum ring radius that leads to this spot-splitting is r = 0.80. Based 
on the relative locations of the equilibrium state and the spot-splitting threshold on the solid curve for Sk = Sk{r) 
in Fig. |44(d)[ we conclude that a dynamically-triggered spot self-replication instability is not possible. 

Next, suppose that A is increased to ^ = 82, and consider an initial ring radius of r = 0.62 at time t = 0. Then, 
from the heavy solid curve in Fig. 44(c) we calculate the initial source strengths as Sc ~ 4.03 on the lower branch 
of the Sc versus r curve, and correspondingly Sk = 3.38 at t = on the upper branch of the Sk versus r curve in 
Fig. |44(d)[ Therefore, we predict that the initial pattern is stable to spot-splitting. However, for t > 0, we predict 
from Fig. |44(d)] that the center-spot will eventually split at the ring radius r = 0.68 when Sk crosses above S2- This 
occurs at a time before the ring radius has approached its equilibrium state at j'ei = 0.71. Therefore, our asymptotic 
theory predicts that a dynamical-triggered spot-splitting instability will occur for this parameter set. This behavior 
is confirmed by the full numerical results shown in Fig. l46l computed from the GS model 



Figure 46. Experiment 1.5: Fix A — 82, D — 3.912, and e = 0.02. The full numerical results computed from the GS 
model Ijl.ip are shown at different instants in time. The initial pattern has nine spots equally spaced on a ring of initial 
radius r = 0.62 together with a center-spot. From left to right the plots of v are shown at times t — 60, 93, 144, 214. 



8 Discussion 

For the GS model ()1.1|) in a two-dimensional domain in the semi-strong interaction regime D = 0(1), and for 
A = ©(— elne), the slow dynamics and instability mechanisms of quasi-equilibrium multi-spot patterns was studied 
by a hybrid asymptotic-numerical method. By using a singular perturbation approach that accounts for all terms 
in powers of ly = — 1/lne, a multi-spot quasi-equilibrium pattern was constructed by asymptotically matching a 
local approximation of the solution near each spot to a global representation of the solution defined in terms of 
the reduced- wave Green's function. The local problem near each spot, referred to as the core problem, consisted of 
a radially-symmetric BVP system that must be solved numerically. In the global, or outer, representation of the 
solution, each spot at a given instant in time is asymptotically approximated as a Coulomb singularity for u of 
strength Sj at location xj e for j = 1, . . . , fc. In Principal Result 3.1 a DAE system was derived for the slow 
time evolution of the coordinates Xj and Sj for j = 1, . . . ,k, which characterizes the slow dynamics of the fc-spot 
quasi-equilibrium pattern. The method used to construct quasi-cquilibria is an extension of that developed in |57j 
(see also [32], |13) . [49]) to sum logarithmic asymptotic expansions for linear PDF models in perforated 2-D domains. 

The stability of the multi-spot quasi-equilibrium solution to three distinct types of 0(1) time-scale instabilities was 
studied. From a numerical study of a local eigenvalue problem near each spot, an explicit criterion for the initiation 
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of a peanut-splitting instability was obtained, and it was verified numerically that this instability triggers a nonlinear 
spot self-replication event. In addition, a singular perturbation method was used to formulate a global eigenvalue 
problem characterizing the onset of either oscillatory or competition instabilities for the multi-spot pattern. This 
global eigenvalue problem was studied for special multi-spot patterns for which a certain Green's matrix is circulant. 

Our hybrid asymptotic-numerical framework for spot dynamics and stability was implemented numerically for some 
spot patterns in M'^ and in the unit disk and square. The results from our theory were shown to very compare favorably 
with those obtained from full numerical simulations of the GS model (II. ip . One key finding is that dynamically- 
triggered instabilities can occur for the GS model in a wide parameter range. The qualitative differences in multi-spot 
solution behavior between the GS model (|l.ip and the Schnakenburg model of |38) were highlighted. 

The numerics required to implement the hybrid asymptotic-numerical theory is rather simple, provided that the 
reduced-wave Green's function for the domain can be readily determined. For the domains considered in this paper, 
this Green's function is available analytically, whereas for other domains it can be computed rapidly by using fast 
multipole methods for linear PDE's (cf. |12j ). When this Green's function is available, the implementation of the 
hybrid asymptotic-numerical approach only requires the numerical solution of a coupled DAE system for the spot 
dynamics and a BVP to compute competition and oscillatory instability thresholds. The DAE system for Xj and Sj 
involves two functions (i.e. xi^j) ^'^d ^{Sj)) that are determined in terms of the solution to the core problem. These 
functions are readily pre-computed by solving some simple EVP's. Similarly, from a BVP eigenvalue problem, the 
condition for the onset of self-replication is readily pre-computed in terms of a threshold value of the source strength 
S,. 

In this way, a similar hybrid approach should be readily applicable to study spot dynamics and spot stability for 
related two-component RD models of the general form 

Vt = e'^Av-v + f{u)vP , TUt = DAu-{u-ua) + e^'^g{u)v"' , x e f7 , (8.1a) 

for p > 0, m > 0, and for certain classes of f{u) and g(u). The key conditions that are needed for our approach are 
that there is a core problem near each spot for which u — >■ and u has a logarithmic growth at infinity in terms of the 
stretched inner variable, and that the the spots are coupled together via a quasi-static linear elliptic concentration 
field for u. Such a generalization of our overall approach for the GS model is undertaken in |10| for the simple case 
of a two-spot evolution in R^. The work in [lOj is the two-dimensional counterpart of the study of [19J of two-spike 
solutions for a general class of RD systems on the infinite line. 

We conclude this paper by discussing a few open problems for the GS model related to spot dynamics, equilibria, 
and stability. The first main problem is to develop a weakly nonlinear stability theory for the three different types of 
instabilities. Our numerical results suggest that both the peanut-splitting instability and the oscillatory instability 
are subcritical, as they lead to spot self-replication and spot-annihilation, respectively. In particular, with regards to 
self-replication, it would be interesting to perform a weakly nonlinear analysis of the time-dependent core problem 
(|4.6p to show that the initial peanut-splitting instability does not saturate in the weakly nonlinear regime, but instead 
triggers a nonlinear spot self-replication event. 

For the equilibrium problem, a main open problem is to use the spot equilibrium condition (j3.8p to calculate 
detailed bifurcation diagrams in terms of A and D of steady-state fc-spot patterns in arbitrary domains. One goal 
would be to determine whether asymmetric equilibria, characterized by equilibrium spot patterns with different 
source strengths, are possible and to then determine the multiplicity of equilibrium solutions. For these asymmetric 
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multi-spot equilibria, it should be possible to calculate oscillatory and competition instability thresholds from our 
global eigenvalue problem. The analysis in this paper has largely been restricted to the case where the Green's matrix 
is circulant, which allows for spot patterns with a common source strength. 

Another open problem is to study weak translational instabilities associated with multi-spot equilibrium patterns. 
Such small eigenvalue instabilities, with growth rate A — C'(e^), are characterized by instabilities of stationary points 
of the DAE system of Principal Result 3.1 governing spot dynamics. Based on the numerical results in Fig. 15 of 
[38j , one simple scenario where such a weak instability occurs is for a one-ring pattern of spots inside a disk. More 
specifically, we expect that an equilibrium pattern consisting of k equally-spaced spots on a ring concentric within 
the unit disk will be weakly unstable with growth rate A — ©(e^) when k exceeds some threshold. Similar thresholds 
for weak instabilities have been studied in |27) for a ring of spots in the context of a simple interacting particle 
system model, and in [7] and [8] for a ring of Eulerian fluid point vortices on the surface of a sphere. 

Finally, it would be very interesting to exhibit a parameter regime for the GS model where the quasi-equilibrium 
spot pattern undergoes repeated episodes of either spot self-replication, when a spot is too far from any of its neigh- 
bors, or spot-oscillation leading to annihilation, which is triggered when two slowly-drifting spots become too closely 
spaced. For D small, but with D >> 0{e'^), such a replication-annihilation "attractor" should be realizable by 
considering the instability thresholds for two-spot interactions in R^. In this paper, a key first step in the construc- 
tion of this "attractor" has been undertaken by showing that the GS model can robustly support the existence of 
dynamically-triggered instabilities for a collection of spots that are initially stable at t — 0. 
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Appendix A The Reduced- Wave and Neumann Green's Function 

In this appendix we give detailed analytical results for both the reduced-wave Green's function and the Neumann 
Green's function for the unit disk and for a rectangle. These results are needed in ^and ^in order to numerically 
implement the analytical theory of ^2Hllfor spot dynamics and instabilities. 



A.l Green's Function for a Unit Disk 

Let 51 = {x : |x| < 1}. Upon introducing polar coordinates {xo,yo) — (po cos ^o, Po sin6'o), the reduced-wave Green's 
function of (|1.2p satisfies 

Gpp + -Gp + ^Ge8-^G^^-S{p~po)d{e~eo), 0<p<l, O<0<27r, (A.l) 

with boundary conditions G{p, 9 + 2tt) = G{p, 6), Gp(l, 0) — 0, and where G is well-behaved as p — > 0. 
To determine G, we first extract the 6 dependence by introducing the complex Fourier series 

G{p, e- po, Oo) = ^ ^^(P' ^o)e-"' , G„(p; po, ^o) = e"'G(p, 0; po, 0o) dO . (A.2) 



n— — oo 
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From (|A.1[) we obtain that Gnpp + P^^Gnp — ri^p^^Gn — D^^Gn = -p^^S{p — po)e™^° on < p < 1, with boundary 
conditions G„p(l; po, Oa) — and Gri(0; po, Oq) < oo. The solution to this problem has the form 

Alln 
A2 



Gn{p; Pa,0o) 



p 



< p < Po , 

/ ( P \ ^"(w) (_p_ 



Po< p<l, 



(A.3) 



where /n('') and Kn{r) are modified Bessel functions of order n. In (|A.3p . the coefficients Ai and A2 are obtained 
by making (?„ continuous at po and by requiring that the jump in the derivative of G„ at po is [G„p]p(, = — ° . 
Then, upon using the Wronskian determinant for Kn and /„, we readily determine Ai and A2, and obtain 



Gn{p\Po-,Oo) = 



P> 



5 



0<p< 1 



(A.4) 



where we have define p< = niin(po,p) and p> — max(po,p). Therefore, the Fourier expansion for the reduced- wave 
Green's function satisfying (jA.ip is 



G(p,^;po,0o) = ^ J2 

2,71 ^ — ' 



P> 



K: 
I: 



p> 



In 



P< 



0<p<l. (A.5) 



The infinite series representation (jA.5|) converges very slowly. To improve the convergence properties, we must 
extract from G the free-space Green's function G/, which has the Fourier representation 



G/(p,0;Po,0o) 



1 

2^ 



X - Xo 



2n ^ ''[VD ' 



n— — 00 



p< 



0<p<l, (A.6) 



where |x — xo| = \/ p^ -\- p^ — 2ppo cos(0 — ^o)- Then, we can decompose G in (|A.5|) as 



G{p,0;po,0o)^—K„ 
Zn 



|x - xol 



D 



- T 

27r ^ 



Po 



< p < 1 . (A.7) 



In (|A.7p . the first term arises from the source at (po, ^o), while the infinite sum represents the effects of the boundary 
conditions. Then, we use K^n{z) = Kn{z) and /_„(z) = /„(z) to further simplify (|A.7p to 



G{p,e;po.e^) = —Ko 

ZTT 



R 



27r r, (l\ '° 



M 



PO 



lo 



K' (^) 

- - V cos{n{9 - 9o)) 

^ T' i — 1 



TT 



Po 



< p < 1 . (A., 



Since K!^ {^~7d) ^ ^^"^ ^" ("v^) ^ ^ '^^'^ integer n > 0, then 6m is a small error term bounded by 



\Sm\<Pm^- 



1 



Po 



--M+1 I'n (^) 

In Table lA 11 we give the number M of Fourier terms in (|A.8I) required to calculate the reduced-wave Green's 
function within a tolerance of 10~^. The series for P^/ is found to converge fairly fast, especially when the singular 
point is not close to the boundary of the unit circle and/or D is not too large. Equation (jA.8|) for G with Pm < 10~^ 
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po p M Pm 



0.8 0.8 31 8.1737e-009 

0.8 0.5 16 6.7405e-009 

0.8 0.2 8 9.5562e-009 

0.5 0.5 11 4.3285e-009 



Table Al. The number of Fourier terms needed to determine G to within a tolerance 10 * for D = 1. 



is the one used in the numerical simulations in JJlHTl From (jA.8[) we can readily extract the self-interaction, or regular 
part, R of the reduced- wave Green's function, as written in (|1.2 h\ . The gradients of R and G that are needed in 
Principal Result 3.1 to determine the motion of a collection of spots are determined numerically from (|A.8p . 



A. 2 Neumann Green's Function for a Unit Disk 



For D ^ 0(1), we recall from (I2.12p that the reduced-wave Green's function and its regular part can be approximated 
by the Neumann Green's function G'^^^(x;xo) and its regular part _R(^)(x;xo) satisfying (j2.13p . For the unit disk, 
this Neumann Green's function and its regular part are defined in cartesian coordinates by (see eq. (4.3) of [321), 



GW(x;xo) = -^ f-ln|x-xo|-ln 



27r 



i?W(xo;xo) = ^ ('-ln(l- 



xo| 



x|xo| 



|xo| 



Xq 
|xo| 



2Uxr 



IxoH-T 



(A.9 a) 
(A.9 b) 



The gradients of these quantities can be obtained by a simple calculation as 



VGW(x;xo) = -^ 



(x- Xq) 



|xo| 



Vi?W(xo;xo)- — 



xol 



l-|xo| 



;-xo|^ x|xo|^-xo 

where the overbar denotes complex conjugate. In fJlHZl these results allow us to explicitly determine the dynamics of 
a collection of spots inside a disk when D is large, as characterized in Principal Result 3.1 of fJH 

When f2 is the unit disk, and for two specific sets of source and observation points in f2, in Table 1X2] we give 
numerical results for the two-term approximation (see equation (j2.12l) ') 



G = 



D 

m 



R^ R = 



D 



(A.IO) 



with \Q\ — n, which relates the reduced- wave and Neumann Green's functions when D is large. The results for G^^\ 
R'^^\ G, and R, were computed from (jA.9p and (jA.Sp . The results in Table 1X2] are believed to be correct to the 
number of digits shown, and confirm that G ~ G and R ~ R when D is large. 



A. 3 Green's Function for a Rectangle 

Next, we calculate the reduced-wave Green's function G{x,y;xo,yo) in the rectangular domain = [0,L] x [0,d]. 
Since the free-space Green's function G/ is given by G/(x; xq) = {2n)~^Ko ( |x — xo|/a/D ) , we can use the method 
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{p,po,9- 9o) 


= (0.3,0.5,0.0) 


{P,PO 


e-eo) = 


(0.8,0.2, 


7V/6) 


D 


G 


G 


R 


R 


G 


G 


R 


R 


0.1 


0.12182 


0.22154 


-.15651 


-.00196 


0.02473 


0.06198 


-.16303 


-.07467 


1.0 


0.48316 


0.50802 


0.24622 


0.28452 


0.33490 


0.34846 


0.19617 


0.21181 


4.0 


1.45566 


1.46295 


1.22822 


1.23945 


1.29922 


1.30339 


1.16243 


1.16674 


10.0 


3.36978 


3.37281 


3.14465 


3.14931 


3.21150 


3.21325 


3.07484 


3.07660 



Table A 2. Comparison of results for the reduced-wave Green's function G and its regular part R for the unit disk 
with the two-term large D approximation G and R, as given in (jA.10[) . The comparison is made at two sets of source 
and observation points in fl. For D large, we note that G ~ G and R fa R, as expected. 



of images to write the solution to (|1.2[) as 

G(x;xo) = 



^ OO OO 4: 

n— — OO m— — OO l—l 



Ko 



D 



^mn = {xq + 2nL, yo + 2md) , x^^^^ = (-xo + 2nL, yo + 2md) , 
^mn = i^o + 2nL, -yo -f 2md) , x^^l = {-xq + 2nL, -yo 2md) . 



(A.ll a) 

(A.116) 
(A.llc) 



Since KQ{r) decays exponentially as r — >■ oo, then the infinite series in (|A.11|) converges fairly rapidly when D = 0(1). 
The regular part of the reduced- wave Green's function together with the gradients of G and R, as required in Principal 
Result 3.1, are readily evaluated numerically from (jA.ll[) . 



A. 4 Neumann Green's Function for a Rectangle 

For D 3> 0{1) is large, we can approximate G and its regular part R by the corresponding Neumann Green's function 
G^^^ and its regular part i?*^^) by the two-term expansion in p.l2p . This latter Green's function, calculated by an 
analytical Ewald-type summation method, was given analytically in formula (4.13) of [38] as 

G(~) (x; xo) = In |x - xol + (x; xo) , (A.12 a) 

where i?*^^^(x;xQ) is given explicitly by 

OO 

i?W(x;xo) = E - - - - ^"C+,+ 111 - g"C+,-l|i - g"C-,+l|i - gX-.-l) 

n=0 

-i'"(n^)-i,|'""-'"=-i' 

Here points in the rectangle are written as complex numbers. In (|A.12 b\ . q = g^^^'^/'*, while the eight complex- valued 
constants z± ± and C±.± a-re given by 

z+,± EE exp{fi{-\x + xo\ + i{y ± yo))/2) , z_,± = exp{fj,{-\x - xo\ -I- i{y ± yo))/2) , (A.12 c) 

C+,± = exp{fi{\x + xo\-2L + i{y± yo))/2) , C-,± = exp(/i(|a; - xo\ - 2L + i{y ± yo))/2) . (A.12 d) 

Here fj, is defined by fj, — 27r/d. The self-interaction term i?(^)(xo; xq), and the gradients of i?'^) and G'-^' can 
be calculated readily from (jA.12p . In particular, for the unit square [0, 1] x [0, 1] with a singularity at the middle 



L 1 max(x,Xo) 1 
'd ~3 L ^ "2 V L2 
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xq = (0.5,0.5) the self-interaction term R^^^q = i?^^^(xo;xo) calculated from (jA.12[) is simply 



4^ - -- E In (1 - 9") + 1^ - ln(2^) , 9 ^ . (A.13) 
TT ^ — ' iz zvr 

n=l 



Appendix B Survey of NLEP Theory: Comparison of Quasi-Equilibria and Stability 

In |60| . the existence and stability of a one-spot pattern to the GS model on the infinite plane M? was studied. For 
this infinite-plane problem, we can set 13 = 1 in (jl.l|) . bi the inner region near the spot, the leading-order analysis 
of [60j in terms of = — 1/ Ine and for A = O (e[— Ine]^/^), showed that u ^ uq, where uq is locally constant near 
the spot. Moreover, v ^ vq — w{p)/{Auq), where w{p) is the radially symmetric ground-state solution of 

ApW - w + = , 0<p<oo; w(0) > , w'(0) = 0; w ^ as p ^ oo , (B.l) 

where ApW = w" p^^w' . For e — >■ 0, it was shown in ^60] that Uq satisfies the quadratic equation 

L r°° 1 

l-uo , L=—- w^pdp, v = -—. (B.2) 

uo V A'^ Jo In £ 

Therefore, the existence condition for a one-spot equilibrium solution is that 

L<- ^ A>Af^ = 2e^-^, bo = J w^{p)pdp. (B.3) 

Define Lq and 7 by Lq = lime_j.o L ~ 0(1) and 7 = Inr = 0(1), so that r = 0{e~'*), and assume that < 7 < 2. 
Then, the following radially symmetric nonlocal eigenvalue problem (NLEP) was derived in equation (5.5) of 



A.^o - ^0 + 2^^o - '^l' V^'^ ^ f-"^t" - , 0<p<oo, (B.4) 
2mo + (1 - uo)(2 - 7) Jo w^pdp 

with V-'o as p — > 00. An analysis of this NLEP in [60] led to Theorem 2.2 of [60] , which we summarize as follows: 

(1) There exists a saddle node bifurcation at Lq — such that there are two equilibrium solutions Uq given by 
Uq — (1 ± Vl ~ 4Lo) /2 for Lq < and no equilibrium solutions when Lq > 

(2) Assume that < 7 < 2 and Lq < -j. Then, the solution branch {uq,Vq) is linearly unstable. 

(3) Assume that < 7 < 2. Then, the other solution branch {uq,Vq) is linearly unstable if 

-1/2 



1 



1 



7 



■7 



A < Agjj} = ^fw 



1 



7 



7 



(B.5) 



(4) Assume that 7 = and Lq < Then, {uq , Wq ) is linearly stable. 

(5) If < 7 < 2, the stability of {uq,Vq) is unknown for A > Agw 

We remark that the saddle-node bifurcation value of A, obtained from (jB.2[) . has the scaling A — O (£[— Ine]^/^) as 
£ — >■ 0. The stability results listed above from [60 also relate only to this range in A. However, the quasi-equilibrium 
spot solution given in Principal Result 2.1 of §2, and the spot self-replication threshold of Principal Result 4.1, 
occur for the slightly larger range A ^ 0(e[— Ine]). Consequently, the occurrence of spot self- replication behavior 
for the GS model (jl.ip was not observed in the scaling regime for A considered in }60| . In §5.1l we compare the Hopf 
Bifurcation threshold predicted by (jB.5l) with that obtained from our global eigenvalue problem of ; 
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B.l NLEP Theory for fc-Spot Patterns 

In [6 2) a related theoretical approach, based on the rigorous analysis of other NLEP's, was used to obtain a leading- 
order asymptotic theory for the existence and stability of fc-spot patterns to the GS model in a bounded 
two-dimensional domain. In the inner region near the spot centered at Xj, the leading-order- in-i/ scaling in j62] 
showed that u ^ Uj, where Uj is a constant, and that v satisfies v ^ vj = -^w{p)^ where 'w{p) is the radially 
symmetric solution of (jB.ip . Moreover, the Uj for j = 1, . . . , fc satisfy the nonlinear algebraic system (cf. |62| ) 

(B.6) 

where Ls and rjs are defined in terms of the area jfij of the domain by 

2e^TTho _^ M , /"^ 2 . 1 



For the case where Uj — uq for j = 1, . . . , fc, then (|B.6[) is a quadratic equation for uj, which leads to the following 
three main parameter regimes of | 



4fcio < 1, forrye^O ^ D:$>0{v-^), 
AijeLe < 1, for T^e ^ oo ^ £» = 
4(fc-|-r/o)ivo < 1, forrye^?7o ■i^ D = 0{i^-^) . 



4?7eie < 1, for T^e ^ oo ^D = 0{1), Lq ^ limLe ^ 0(1) , 770 = hni ?]£ = 



(B.8) 

From (|B.7|) . the second line in ()B.8P is for A = O Ine]^/^) with D = 0(1), while the third line in (|B.8|) is for 
A ~ O (e) with D = 0{v~^). v — — 1/lne. Thus, the fc-spot quasi-equilibrium patterns of }62j do not correspond to 
the range A = O {e[— Ine]) and D — 0(1), considered in ^JSJ where spot-splitting can occur. 

As shown in (62j . the leading-order quasi-equilibrium solution satisfies Uj ^ Uq in the inner region. The global 
representation for v is v Vq . Here Uq and satisfy 

■i (1 ± VI - 4fcio) , for i?»0(z/-i), ^ 
u^^\ i (1 ± VI - 4^VeLs) , fori? = 0(l), ^± ^ ^ _i_ (^-ijx - x, |) . (B.9) 

i (1 ± v/1 - 4(fc + r;o)Lo) , for i? = 0(i^-i), 

The stability analysis in |62j was largely based on the derivation and analysis of the radially symmetric NLEP 



Ap^/) - -0 + 2wV' - 7 w — 7; ^ 1 0<p<oo; ip ^ , p-^oo. (B.lOa) 

Jo w^pdp 

In terms of uq — Uq and r — 0(1), in (|B.10 there are two choices for 7 given by 



2r]£L£{l + TX) +kL£ 27]£L 



e 



{ul + L£T]£)(1 + t\) + kLe ' ul + Lsrje 



(B.10&) 



Assuming that the third condition in (|B.8[) holds. Theorem 2.3 of |62j gives a stability result for fc-spot equilibrium 
solution based on an analysis of the NLEP (|B.10[) . together with certain properties of J^(xi, . . . , x^) defined by 



k k 

^(xi,...,Xfc)=-^^(C?),,^. . (B.U) 

i=i j=i 



Here Q is the Green's matrix of ()2.9p defined in terms of the reduced- wave Green's function. From [62], an equilibrium 
spot pattern must be at a critical point of T. We denote by "Ho the Hessian of J- at this critical point. 
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Suppose that A{r]s + k)Ls < 1/4. Then, Theorem 2.3 of |62| proves that the large solutions are all hnearly 

unstable. For the small solutions (u^, w"), this theorem provides the following results for three ranges of rye: 

(1) Assume that rje 0, so that D > ©(i^^^). Then, 

• If fc = 1, and all the eigenvalues of the Hessian "Hq are negative, then there exists a unique ti > 0, with 
Ti — 0(1), such that for r < ri, {u^,v^) is linearly stable, while for r > ti, it is linearly unstable. 

• If fc > 1, is linearly unstable for any r > 0. 

• If the Hessian Ho has a strictly positive eigenvalue, then {u" ,v~) is linearly unstable for any r > 0. 

(2) Assume that rie -> oo, so that D < 0{iy~'^). Then, 

• If all the eigenvalues of the Hessian Hq are negative, then , v") is linearly stable for any r > 0. 

• If the Hessian Hq has a strictly positive eigenvalue, {u~ ,v~) is linearly unstable for any r > 0. 

(3) Assume that 1]^ t/q, so that D = 0{h'~^). The following results hold: 

• If Lo < (2rjo+ky^ ' ^'^'^ eigenvalues of Hq are negative, then {u^,v^) is linearly stable for r = 0(1) 
sufficiently small or r = 0(1) sufficiently large. 

• If fc = 1, Lq > (2i)o+i)^ ' ^^'^ ^^"^ eigenvalues of Ho are negative, then there exists T2 > 0, T3 > 0, such 
that for T < T2, {u^,v^) is linearly stable, while for r > T3, it is linearly unstable. 

• If fc > 1 and Lq > (2?)o+fc)^ ' then (w~,w~) is linearly unstable for any r > 0. 

• If the Hessian Hq has a strictly positive eigenvalue, {u~ , v") is linearly unstable for any r > 0. 

We now relate these results of ^62j with those obtained in ^IHH Firstly, the condition of |62j that an equilibrium 
/c-spot configuration Xi , . . . , x^: must be at a critical point of T in ()B.11|) is equivalent to the equilibrium result in 
p.Sp of Principal Result 3.1 provided that we replace the source strengths Sje in p.Sp with their leading-order- 
in- 1/ asymptotically common value Sc ^ A, as obtained from (|2.1ip . The stability condition of [62 expressed in 
terms of the sign of the eigenvalues of the Hessian Ho of J- is equivalent to the statement that the equilibrium spot 
configuration Xi, . . . , x^ is stable with respect to the slow motion ODE dynamics of Principal Result 3.1 of fj3] when 
we use the leading-order approximation Sj ^ Sc for j = 1, . . . , fc. Therefore, the condition on the Hessian Ho in [62j 
relates to the small eigenvalues of order O(e^) in the linearization. 

In addition, the second result of Theorem 2.3 of [62j for r/s — > 00, as written above, includes the range D = 0(1). 
The condition 4(77^ -I- k)Ls < 1/4, together with (|B.7|) relating Le to A, shows that this second result of [6 2) holds 
when A = O (^e[— Ine]^^^^ . Therefore, for this range of A and D, the leading order NLEP analysis in |62| predicts 
that there is an equilibrium solution branch whose stability depends only on the eigenvalues of the Hessian matrix 
Ho- This result is very similar to that obtained from our global eigenvalue problem in Principal Result 4.2 of 



(see Fig. 9(a) ) for the nearby parameter range A = 0(— e Ine). There, we showed to leading-order in ly that a /c-spot 
quasi-equilibrium solution is stable to competition or oscillatory instabilities when D = 0(1) and A = O (— elne). 

With the exception of the condition on the Hessian matrix of J^, these stability results of |62j depend only on the 
number k of spots, and not on their spatial locations within il. This is qualitatively very different to that obtained 
in Sj5]-Sj5] from our globally coupled eigenvalue problem, which accounts for all terms in powers of v. 



(B.12) 
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B.2 Comparison of the Quasi-Equilibrium Solutions 

In this subsection wc show that our asymptotic construction of quasi-cquiUbrium fc-spot solutions, as given in Principal 
Result 2.1 of 351 can be reduced to that given in [62] when A = 0{v^/'^) in ^11^. Note that when A = 0(i^^/^), then 
from (|2.5p . O {^e[— Ine]^^^^ , which is the key parameter regime of j62j . 

For A — 0{v^/'^), ()2.7p together with the core problem ()2.3|) . suggest that we expand Uj, Vj, Sj, and x, in ()2.3|) as 

X - J^"'^'(X0j + J^Xij + ■••), v-^'^iUa, +vUi, + ■■■), 

S.^v^l^ (5oj + vSx, + •••), V-j^v^l^ ( + t^Vij + • • • ) • 
Upon substituting (jB.12p into p.3p . and collecting powers of v, we obtain the leading-order problem 

AC/oj=0, C/oj^Xoj asp^oo; AVbj - Fqj + C/oj^oj = , Foj ^ as p ^ oo . (B.13) 

The solution to (IB.ISP is C/qj — Xoj: and Vqj = wjxoj^ where wi^p) satisfies (IB.1[) . At next order, we get 

MJij = UojVQj , Uij ^ Sqj In p + Xij" as p oo , (B.14 a) 

AVij - Vij + 2UojVojVij ^ -UijV^^ , Vij^O as p ^ oo . (B.14&) 

At one higher order, we find that U2j satisfies 

AU2j ^ 2UojVojVij + V^jUij , U2j Sij\np + X23 as p ^ oo . (B.15) 

Upon applying the divergence theorem to (|B.14 a[) . and using Uoj — Xoj and Voj = w/xoji we obtain 

/■CO ^ /■CO 

Soj ^ / UojVi^j p dp = , bo= / w'^pdp, (B.16) 

JQ XOj Jq 

which determines Sqj in terms of xoj ■ We then decompose the solution to (jB.14[) in the form 



^ (xojXij + Uij) , "^ij = -4- ( - XojXij w + Vij) 



XO] ^ ^ Xoj 

where Uij and Vij are the unique solutions on < p < oo of 

ApUij=w'^, ApVij - Vij -\- 2wVij ^ --Uij w'^ ; Vij ^ Uij - bo\n p , as p oo . (B.17) 
Similarly, by applying the divergence theorem to (jB.15[) . we calculate Sij as 

5i, = + , bi= f {w^Uij + 2wVij)pdp. (B.18) 

Xoj Xoj J 

The BVP solver COLSYS (cf. [2]) is used to numerically compute the ground-state solution w{p) together with the 
solutions Uij and Vij of (|B.17|) . Then, from a simple numerical quadrature, we estimate bo « 4.9347 and bi « 0.8706. 
Finally, upon substituting (jB.16P and (jB.lSP into the nonlinear algebraic system p.7p of Sj2l we obtain that 



A 1/2 1 I ^0 \ , 3/2 r ^1 I / 1 ^0 

V Xoj/ Lxoj Y Xoj 



Xij + 27r 




(B.19) 



For a prescribed A with A ~ Anv^l'^ + Aiv^^^ + • • • , we could then calculate xoj and xij from (jB.lQp . Since 
A = 0(y^l'^) corresponds to A = O [— Ine]^^^^, we conclude that it is in this parameter range that the coupled 
core problem (|2.3p reduces to the scalar ground-state problem for w, as was exploited in [62] . 

Finally, we show that (IB.19P reduces to the expression (jB.6|) used in |62) . To see this, we note from (|2.ip that 
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A = A^/Dvje and = eUj|{A^JD) « ev~^l'^xni,l{A^)- Moreover, for D > 0(1), we use = |n[ + ^(1)' and 



= -p- + 0(1) for i 7^ j (see equation (|2?T2l) of Then, from (|B.19p . we obtain 



A^JDv As/Dv 



A-JD Ui 



2tivD 



boe 



AvDui A^/Du, 



Upon using (|B.7p for Ls and 775, (jB.20p reduces to the following result, in agreement with (jB.6P of }62j : 



1 — Ui 



boe' 



A^Di 



+ 



27r 6o£^ 1 Ls-qe 



M A^ 



E- 



(B.20) 



(B.21) 



2 = 1 



B.3 Comparison of the Global Eigenvalue Problem with the NLEP of 

Next, we show for A = OijAl'^) and D = 0[v^^) that our eigenvalue problem, consisting of (|4.8p coupled to (|4.13p . 



can be reduced to the NLEP (|B.10p derived in |62] . To show this, we expand $j, Nj, Bj, and Cj, in ()4.8p as 

Bj = Bjo + vBji H , = u{Cjo + vCji H ) , N.j^N.jo + vNji H , + 

Upon substituting these expansions, together with Uj ~ v^^/'^xoj a-nd Vj ~ v^l'^w jxQji into (|4.8p for A'j, we conclude 

that A^jo = ^jo, and that A'^i satisfies 

9 

ApA'ji = —Bjo + 2w'ipj , < p < 00 ; Ar,i ^ Qo In p + Bji as p ^ cxi . (B.22) 
Xoj 

Upon applying the divergence theorem to (jB.22[) . we calculate Cjo as 



bo f°° 

X01 Jo Jo 



w pdp . 



(B.23) 



Xoj Jo Jo 

Then, we let D = Dq/v ^ 1 where Dq = 0(1), and we impose that xoj = Xo for j = 1, . . . , fc, so that to leading- 
order each local spot solution is the same. For Z) 3> 1, the A-dependcnt Green's function G\ij and its regular part, 
R\ jj , satisfying ()4.1ip , have the leading-order behavior 

D , , „ D 



G 



Xij 



+ 0(1), i^j; Rxj„ 



|r!|(l + TA) • ^ ' |r!|(l + rA) 

Upon using (|B.24p . together with Bj ^ Bjq and Cj ^ i^Cjo, we obtain that (|4.13l) . with D = Dq/i^, reduces to 

2ttD, 



(B.24) 



{l + rXm ^ 

Upon substituting (|B.23P for Cjo into this equation, we obtain that for j = 1, . . . , satisfies 



Xo 



2^Da bo 



1^ pOO 



AnDo 



{i + TX)\n\ ^ 



^ /"OO 

V / wij.pdp. (B.25) 

2 = 1-^0 



Next, we define the vectors po and by po = (Bio, • • ■ , Bko)^ and ^ = (/J^ pdp,. . . , ui^pk pdp) ■ Then, 
(|B.25p can be written in matrix form in terms of e = (1, . . . , 1)-^ as 



1 + ^1/ 



2t:D, 



"0 T 

ee 



Po = - 2/ 



4ttD, 







(l + rA)|r! 



■ee^ I ^ . 



xi) {i+T\mxi 

The matrix multiplying po is an invertible rank-one perturbation of the identity matrix, and hence 

-1 



V : 



1+^1/ 



27rD, 



oo T 
ee 



Xl) [l + rXmxl 



21 



(l + rA)|r!| 



T 

ee 



(B.26) 
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We then re- write some quantities in the symmetric matrix D in terms of Lg and 77^, as defined in (|B.7[) . to get 



Finally, to obtain an NLEP we substitute Uj ~ v ^^^Xo, Vj ^^^^^w/xo, i^^j, and Njq ^ Bjq, into (I4.8P for 

^j, which leads to the following radially symmetric vector NLEP for 1/) = {ipi, . . . ,ipk)'^: 



where ri and r2 are the eigenvalues of T) given in (|B.27p . This NLEP is identical to the NLEP (jB.lOp studied 
rigorously in [62] . Therefore, we conclude that our global eigenvalue problem in ii4.2l reduces to leading-order in v to 
the NLEP problem of [62] when A = o(e[-~ Ine]^/^) and D = 0{v-^). 



Appendix C Numerics for the Global Eigenvalue Problem: Competition and Oscillatory Instabilities 

In [J5]-[J6] stability thresholds are computed for spot configurations Xi, . . . , x^ for which the Green's matrix Q and the 
A— dependent Green's matrix Qx are circulant. Therefore, we must solve the coupled eigenvalue problem (|4.17l) and 
(|4.14p . where the common spot profile and source strength Sc satisfy p.3p and (I2.22p . respectively. In solving (|2.3p 
and (|4.14p . we have chosen the finite interval [0, L] with L = 15 to adequately approximate the infinite domain. 

One numerical approach is to calculate the eigenvalue A by a fixed-point iterative type method for a given set of 
parameters A, D, r, and £, and for a given configuration of spots. The outline of this method is as follows: 

(1) Given T,e,A, and D, calculate Sc, x{Sc), and the common spot profile Uc,Vc from p.3p and (|2.22p . 

(2) For the rfi^ iteration, starting from a known eigenvalue A", solve (I4.14p for the common value Be- 

(3) Next, with Be known, compute a new approximation A"^^ to A by numerically solving ()4.17p using Newton's 
method. This requires the evaluation of the matrix eigenvalues of the A-dependent Green's matrix in (j4.20l) . 

(4) Repeat step 2 using the updated value A"+^ until the algorithm converges to some A, with an error that is less 
than a given tolerance. Then, A is the eigenvalue of the globally coupled eigenvalue problem. 

This fixed-point iteration approach for A converges rather slowly, and its success relies on the initial guess. In 
addition, since we need to study the stability for a range of values of the parameter r, the computational effort 
required with this simple scheme is rather large. In ^SHSl this scheme is used to plot the path of certain eigenvalues 
in the complex plane as a parameter is varied. 

Our primary numerical method used in ^JSHS] is based on the assumption that an instability occurs at some critical 
values of the parameters. For instance, to compute the Hopf bifurcation threshold th, we assume that there is a 
complex conjugate pair of pure imaginary eigenvalues when t = th for which A^ = Re(A) = and A^ = Im(A) ^ 0. 
Instead of solving for A for each given r, we fix A^ = 0, and solve for the thresholds th and A^, associated with a 
specific eigenvalue uj\j and eigenvector oi Q\. A rough outline of the algorithm that we use is as follows: 




A„V — "0 + 2^0 H o — ^735 = A'0, 0<p<oo; -0—^0, asp— >-cx). 

Xo Jo pw^dp 

By diagonalizing this vector NLEP by using the eigenpairs of 2?, we obtain the scalar NLEP 





< p < 00 ; 



(B.28) 
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(1) Given e, A, and D, calculate 5'c, x(S'c), and the common spot profile Uc^ Vc from (|2.3p and (|2.22p . 

(2) Fix Xr = 0. In the rfi^ iteration, starting from the current approximation t^"\ solve the BVP (|4.14p with 
boundary conditions N^{L) = 1/L and $^(L) = using COLSYS (cf. [2j). This yields Be as Be = Nc{L)-\nL, 
and its derivative dx-Bc is calculated numerically by varying Xi using a centered difference scheme. 

(3) Calculate Q\ and evaluate its 

jth ei genvalue ujj{TXi). The partial derivatives dr(^j and d\^ujj are computed 
numerically by a centered difference scheme. 

(4) Calculate the residuals Re(/j"'') and Im(/j"'') in (|4.18p . and then calculate the Jacobian 



J ■ 



9.Re(/j"^) aA.Re(/j"M ^ / a,Re(4) dx^Ke{B,) + 2T:dx,Ke{u:j) 
ajm(/j")) dxMiff) ) \drlHBc) dxM{Bc) + 2^dxM{^o) 



(5) Use Newton's method to update (n+i) ] — \ (n) ^\ . Then go to step 2 and iterate 

further until reaching a specified tolerance. 

To compute competition instability thresholds, we can use either the simpler formulation (|4.24[) . or else employ a 
similar approach to that outlined above for computing oscillatory instability thresholds. The competition instability 
threshold occurs at A = when the largest eigenvalue first enters the right half plane along the Im(A) ~ axis. 
For this case, the A— dependent matrix satisfies Qx=o = and both Be and ujj are real-valued. We fix all the other 
parameters, and choose either A, D or the spot location as the bifurcation parameter. For instance, if we treat A as 
the bifurcation parameter, then the algorithm to compute the competition stability threshold is as follows: 

(1) Set A = and fix all of the parameters except for A. Calculate Q and its eigenvalue ujj. 

(2) In the n*^ iteration, starting from initial guess for A, we calculate Sc and Uc, Vc from (|2.22p and (|2.3p . 

(3) Solve the BVP (jiH)) with A ^ and N'^{L) = 1/L and ^'^{L) = 0. Compute B^ as Be = Nc{L) - InL. 

(4) Substitute Be and ujj in (|4.17p . and calculate the real- valued residual /j"''. 

(5) Compute numerically by a centered difference scheme. Also compute —^j- = ~ "ff^ 

(6) In gUl), use Newton's method to update = - (^9/j"V<9^) /]"''■ Then, go to step 2 and iterate 
further until reaching a specified tolerance. 

In computing the thresholds for an oscillatory instability for the infinite-domain problem in Sj5l and for the unit 
square and disk in f|6l we must evaluate the A-dependent Green's matrix Qx- This is done by replacing D with 
D/{1 + rA) in the explicit formulae of Appendix A. However, since A is in general complex, in order to use the 
results in Appendix A we must be able to evaluate the required modified Bessel functions Kn{z) and In{z) of a 
complex argument. This was done using the special function software of |64]. 

In ^5]-^ the numerical approach outlined in this appendix was used to compute oscillatory and competition 
instability thresholds for various spot configurations in the infinite plane as well as the unit square and disk. 
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